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

Merge remote-tracking branch 'upstream/develop' into diff_tally3

Browse source
This commit is contained in:
Sterling Harper 2016-01-18 16:25:38 -05:00
commit e2c2f8c939
157 changed files with 5941 additions and 3256 deletions
Download patch file Download diff file Expand all files Collapse all files

6
.gitignore vendored
View file

@ -15,8 +15,9 @@ openmc.egg-info/
# Inputs generated from Python API
examples/python/**/*.xml
# emacs backups
# emacs and vim backups
*~
*.swp
# OpenMC statepoints
*.binary
@ -24,6 +25,7 @@ examples/python/**/*.xml
# Documentation builds
docs/build
docs/source/_images/*.pdf
docs/source/_images/*.aux
# Source build
build
@ -73,4 +75,4 @@ docs/source/pythonapi/examples/*.xls
docs/source/pythonapi/examples/mgxs
docs/source/pythonapi/examples/tracks
docs/source/pythonapi/examples/fission-rates
docs/source/pythonapi/examples/plots
docs/source/pythonapi/examples/plots

2
CMakeLists.txt
View file

@ -263,7 +263,7 @@ endif()
# set compile flags. Note that this sets the COMPILE_OPTIONS property (also
# available only in 2.8.12+) rather than the COMPILE_FLAGS property, which is
# deprecated. The former can handle lists whereas the latter cannot.
if(CMAKE_VERSION VERSION_LESS 4.8.12)
if(CMAKE_VERSION VERSION_LESS 2.8.12)
string(REPLACE ";" " " f90flags "${f90flags}")
set_property(TARGET ${program} PROPERTY COMPILE_FLAGS "${f90flags}")
else()

9
docs/Makefile
View file

@ -17,6 +17,8 @@ ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) sou
SVG2PDF = inkscape
PDFS = $(patsubst %.svg,%.pdf,$(wildcard $(IMAGEDIR)/*.svg))
# Tikz to PNG conversion
PNGS = $(patsubst %.tex,%.png,$(wildcard $(IMAGEDIR)/*.tex))
.PHONY: help images clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest
@ -43,8 +45,13 @@ help:
%.pdf: %.svg
$(SVG2PDF) -f $< -A $@
%.png: %.tex
pdflatex --interaction=nonstopmode --output-directory=$(IMAGEDIR) $<
pdftoppm -r 120 -singlefile $(patsubst %.tex,%.pdf, $<) $(basename $<)
convert -trim -fuzz 2% -transparent white $(patsubst %.tex,%.ppm,$<) $@
# Rule to build PDFs
images: $(PDFS)
images: $(PDFS) $(PNGS)
clean:
-rm -rf $(BUILDDIR)/*

2
docs/requirements-rtd.txt Normal file
View file

@ -0,0 +1,2 @@
sphinx-numfig
jupyter

BIN
docs/source/_images/cmfd_flow.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 25 KiB

12
docs/source/methods/cmfd_tikz/cmfd_flow.tikz → docs/source/_images/cmfd_flow.tex
View file

@ -1,3 +1,12 @@
\documentclass{standalone}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.11}
\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy}
\pagestyle{empty}
\begin{document}
\begin{tikzpicture}
\matrix[every node/.style={draw, thick, minimum width=3cm, minimum height=1cm, align=center}, column sep=2cm, row sep=1cm] (m) {
\node[draw, fill=red!40] (start) {Batch $i$ \\ tally NDA}; & \\
@ -16,4 +25,5 @@
\draw (modify.north) -- (end.south);
\end{scope}
\end{tikzpicture}
\end{tikzpicture}
\end{document}

BIN
docs/source/_images/meshfig.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 46 KiB

41
docs/source/methods/cmfd_tikz/meshfig.tikz → docs/source/_images/meshfig.tex
View file

@ -1,10 +1,20 @@
\documentclass[tikz]{standalone}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.11}
\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy}
\usepackage{fixltx2e}
\pagestyle{empty}
\begin{document}
% these dimensions are determined in arrow_dimms.ods
\def\scale{1.0}
\def\latWidth{0.2808363589*\scale}
\def\RPVOR{3*\scale}
\def\rectW{0.75*\scale}
\def\RPVIR{2.8694005485*\scale}
@ -20,7 +30,7 @@
\def\bafMIRy{1.5445999739*\scale}
\def\bafMORx{1.8544620609*\scale}
\def\bafMORy{1.573625702*\scale}
\tikzset{Assembly/.style={
inner sep=0pt,
text width=\latWidth in,
@ -29,13 +39,13 @@
align=center
}
}
\def\tkzRPV{(0,0) circle (\RPVIR) (0,0) circle (\RPVOR)}
\def\tkzBarrel{(0,0) circle (\BarrelIR) (0,0) circle (\BarrelOR)}
\def\tkzShields{(0,0) circle (\BarrelOR) (0,0) circle (\ShieldOR)}
\def\tkzBaffCOR{(-\bafCORx, -\bafCORy) rectangle (\bafCORx, \bafCORy)}
\def\tkzBaffCIR{(-\bafCIRx, -\bafCIRy) rectangle (\bafCIRx, \bafCIRy)}
\def\tkzBaffCIR{(-\bafCIRx, -\bafCIRy) rectangle (\bafCIRx, \bafCIRy)}
\def\tkzBaffMOR{(-\bafMORx, -\bafMORy) rectangle (\bafMORx, \bafMORy)}
\def\tkzBaffMIR{(-\bafMIRx, -\bafMIRy) rectangle (\bafMIRx, \bafMIRy) }
\def\tkzBaffleC{ \tkzBaffCIR \tkzBaffCOR }
@ -53,7 +63,7 @@
\begin{tikzpicture}[x=1in,y=1in, xshift=3in]
\scalebox{0.6}{
% draw RPV, barrel, and shield panels
\path[fill=black,even odd rule] \tkzRPV;
\path[fill=black,even odd rule] \tkzBarrel;
\begin{scope}
@ -61,9 +71,9 @@
\path[fill=black,even odd rule] \tkzShields;
\end{scope}
% draw assembly row/column headers
\draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {R} -- ($(-7*\latWidth,4*\latWidth)$);
\draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {P} -- ($(-6*\latWidth,6*\latWidth)$);
\draw[red, thick] ($(-5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {N} -- ($(-5*\latWidth,7*\latWidth)$);
@ -79,7 +89,7 @@
\draw[red, thick] ($(5*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {C} -- ($(5*\latWidth,7*\latWidth)$);
\draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {B} -- ($(6*\latWidth,6*\latWidth)$);
\draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[above, anchor=south] {A} -- ($(7*\latWidth,4*\latWidth)$);
\begin{scope}[rotate=90]
\draw[red, thick] ($(-7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {15} -- ($(-7*\latWidth,4*\latWidth)$);
\draw[red, thick] ($(-6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {14} -- ($(-6*\latWidth,6*\latWidth)$);
@ -97,7 +107,7 @@
\draw[red, thick] ($(6*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {2} -- ($(6*\latWidth,6*\latWidth)$);
\draw[red, thick] ($(7*\latWidth,\RPVOR/\latWidth*\latWidth)$) node[left, anchor=east] {1} -- ($(7*\latWidth,4*\latWidth)$);
\end{scope}
% draw fuel assembly nodes
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,8*\latWidth)$) {};
@ -116,7 +126,7 @@
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,7*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,7*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,7*\latWidth)$) {};
@ -581,7 +591,7 @@
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 6*\latWidth,-7*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 7*\latWidth,-7*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-7*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-8*\latWidth,-8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-7*\latWidth,-8*\latWidth)$) {};
\node [Assembly, fill=\lightgray, opacity=0.3] at ($(-6*\latWidth,-8*\latWidth)$) {};
@ -601,7 +611,7 @@
\node [Assembly, fill=\lightgray, opacity=0.3] at ($( 8*\latWidth,-8*\latWidth)$) {};
% draw baffle north/south
\begin{scope}[even odd rule]
\clip[rotate=90] \tkzBaffMClip;
\path[fill=black] \tkzBaffleC;
@ -611,9 +621,9 @@
\clip \tkzBaffMClip;
\path[fill=black, rotate=90] \tkzBaffleM;
\end{scope}
% draw baffle east/west
\begin{scope}[rotate=90]
\begin{scope}[even odd rule]
\clip[rotate=90] \tkzBaffMClip;
@ -626,3 +636,4 @@
\end{scope}
\end{scope}}
\end{tikzpicture}
\end{document}

BIN
docs/source/_images/openmc200px.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 6.8 KiB

46
docs/source/conf.py
View file

@ -13,6 +13,26 @@
import sys, os
# Determine if we're on Read the Docs server
on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
# On Read the Docs, we need to mock a few third-party modules so we don't get
# ImportErrors when building documentation
try:
from unittest.mock import MagicMock
except ImportError:
from mock import Mock as MagicMock
class Mock(MagicMock):
@classmethod
def __getattr__(cls, name):
return Mock()
MOCK_MODULES = ['numpy', 'h5py', 'pandas', 'opencg']
sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
@ -26,9 +46,8 @@ sys.path.insert(0, os.path.abspath('../..'))
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = ['sphinx.ext.autodoc',
'sphinx.ext.napoleon',
'sphinx.ext.pngmath',
'sphinx.ext.mathjax',
'sphinx.ext.autosummary',
'sphinxcontrib.tikz',
'sphinx_numfig',
'notebook_sphinxext']
@ -105,16 +124,15 @@ pygments_style = 'tango'
# The theme to use for HTML and HTML Help pages. Major themes that come with
# Sphinx are currently 'default' and 'sphinxdoc'.
html_theme = 'haiku'
#html_theme = 'altered_nature'
#html_theme = 'sphinxdoc'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
html_theme_options = {'full_logo': True,
'linkcolor': '#0c3762',
'visitedlinkcolor': '#0c3762'}
if on_rtd:
html_theme = 'default'
html_logo = '_images/openmc200px.png'
else:
html_theme = 'haiku'
html_theme_options = {'full_logo': True,
'linkcolor': '#0c3762',
'visitedlinkcolor': '#0c3762'}
html_logo = '_images/openmc.png'
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = ["_theme"]
@ -126,10 +144,6 @@ html_title = "OpenMC Documentation"
# A shorter title for the navigation bar. Default is the same as html_title.
#html_short_title = None
# The name of an image file (relative to this directory) to place at the top
# of the sidebar.
html_logo = '_images/openmc.png'
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.

11
docs/source/devguide/docbuild.rst
View file

@ -12,15 +12,8 @@ is via pip:
sudo pip install sphinx
Additionally, you will also need two Sphinx extensions for TikZ support and
numbering figures. The sphinxcontrib-tikz_ package should be installed directly
from the git repository as such:
.. code-block:: sh
sudo pip install https://bitbucket.org/philexander/tikz/get/HEAD.tar.gz
The Numfig_ package can be installed directly with pip:
Additionally, you will also need a Sphinx extension for numbering figures. The
Numfig_ package can be installed directly with pip:
.. code-block:: sh

18
docs/source/methods/cmfd.rst
View file

@ -110,10 +110,12 @@ illustrated as a flow chart below. After a batch of neutrons
is simulated, NDA can take place. Each of the steps described above is described
in detail in the following sections.
.. tikz:: Flow chart of NDA process. Note "XS" is used for cross section and
"DC" is used for diffusion coefficient.
:libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy
:include: cmfd_tikz/cmfd_flow.tikz
.. figure:: ../_images/cmfd_flow.png
:align: center
:figclass: align-center
Flow chart of NDA process. Note "XS" is used for cross section and "DC" is
used for diffusion coefficient.
Calculation of Macroscopic Cross Sections
-----------------------------------------
@ -422,9 +424,11 @@ during the MC simulation with incoming and outgoing partial currents. This
allows the user to not have to worry about neutrons producing adequate tallies
in mesh cells far away from the core.
.. tikz:: Diagram of CMFD acceleration mesh
:libs: shapes, snakes, shadows, arrows, calc, decorations.markings, patterns, fit, matrix, spy
:include: cmfd_tikz/meshfig.tikz
.. figure:: ../_images/meshfig.png
:align: center
:figclass: align-center
Diagram of CMFD acceleration mesh
During an MC simulation, CMFD tallies are accumulated. The basic tallies needed
are listed in Table :ref:`tab_tally`. Each tally is performed on a spatial and

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

@ -1626,7 +1626,7 @@ another.
.. _ENDF-6 Format: http://www-nds.iaea.org/ndspub/documents/endf/endf102/endf102.pdf
.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721_3rdmcsampler.pdf
.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf
.. _LA-UR-14-27694: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-14-27694
@ -1636,4 +1636,4 @@ another.
.. _lectures: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-05-4983.pdf
.. _MCNP Manual: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/MCNP5_Manual_Volume_I_LA-UR-03-1987.pdf
.. _MCNP Manual: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-03-1987.pdf

2
docs/source/methods/random_numbers.rst
View file

@ -70,5 +70,5 @@ the idea is to determine the new multiplicative and additive constants in
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
.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf
.. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator

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

@ -507,6 +507,6 @@ improve the estimate of the percentile.
.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution
.. _unpublished rational approximation: http://home.online.no/~pjacklam/notes/invnorm/
.. _unpublished rational approximation: https://web.archive.org/web/20150926021742/http://home.online.no/~pjacklam/notes/invnorm/
.. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf

84
docs/source/pythonapi/examples/mgxs-part-i.ipynb
View file

@ -146,6 +146,8 @@
"\n",
"import openmc\n",
"import openmc.mgxs as mgxs\n",
"from openmc.source import Source\n",
"from openmc.stats import Box\n",
"\n",
"%matplotlib inline"
]
@ -341,7 +343,8 @@
"settings_file.particles = particles\n",
"settings_file.output = {'tallies': True, 'summary': True}\n",
"bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n",
"settings_file.set_source_space('fission', bounds)\n",
"settings_file.source = Source(space=Box(\n",
" bounds[:3], bounds[3:], only_fissionable=True))\n",
"\n",
"# Export to \"settings.xml\"\n",
"settings_file.export_to_xml()"
@ -420,24 +423,22 @@
"data": {
"text/plain": [
"OrderedDict([('flux', Tally\n",
"\tID =\t10000\n",
"\tName =\t\n",
"\tFilters =\t\n",
" \t\tcell\t[1]\n",
" \t\tenergy\t[ 0.00000000e+00 6.25000000e-07 2.00000000e+01]\n",
"\tNuclides =\ttotal \n",
"\tScores =\t['flux']\n",
"\tEstimator =\ttracklength\n",
"), ('absorption', Tally\n",
"\tID =\t10001\n",
"\tName =\t\n",
"\tFilters =\t\n",
" \t\tcell\t[1]\n",
" \t\tenergy\t[ 0.00000000e+00 6.25000000e-07 2.00000000e+01]\n",
"\tNuclides =\ttotal \n",
"\tScores =\t['absorption']\n",
"\tEstimator =\ttracklength\n",
")])"
" \tID =\t10000\n",
" \tName =\t\n",
" \tFilters =\t\n",
" \t\tcell\t[1]\n",
" \t\tenergy\t[ 0.00000000e+00 6.25000000e-07 2.00000000e+01]\n",
" \tNuclides =\ttotal \n",
" \tScores =\t['flux']\n",
" \tEstimator =\ttracklength), ('absorption', Tally\n",
" \tID =\t10001\n",
" \tName =\t\n",
" \tFilters =\t\n",
" \t\tcell\t[1]\n",
" \t\tenergy\t[ 0.00000000e+00 6.25000000e-07 2.00000000e+01]\n",
" \tNuclides =\ttotal \n",
" \tScores =\t['absorption']\n",
" \tEstimator =\ttracklength)])"
]
},
"execution_count": 13,
@ -516,10 +517,9 @@
"\n",
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.0\n",
" Git SHA1: c4b14a5ef87f004528d35cbf33fef3ed15a386ca\n",
" Date/Time: 2015-12-02 09:11:05\n",
" MPI Processes: 1\n",
" Version: 0.7.1\n",
" Git SHA1: ea9fb637f63f9374c7436456141afa850b84acf9\n",
" Date/Time: 2016-01-14 07:16:05\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
@ -604,20 +604,20 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 4.1700E-01 seconds\n",
" Reading cross sections = 8.9000E-02 seconds\n",
" Total time in simulation = 1.4728E+01 seconds\n",
" Time in transport only = 1.4712E+01 seconds\n",
" Time in inactive batches = 1.7890E+00 seconds\n",
" Time in active batches = 1.2939E+01 seconds\n",
" Time synchronizing fission bank = 5.0000E-03 seconds\n",
" Sampling source sites = 3.0000E-03 seconds\n",
" SEND/RECV source sites = 2.0000E-03 seconds\n",
" Time accumulating tallies = 1.0000E-03 seconds\n",
" Total time for initialization = 1.1720E+00 seconds\n",
" Reading cross sections = 9.0300E-01 seconds\n",
" Total time in simulation = 1.7319E+01 seconds\n",
" Time in transport only = 1.7310E+01 seconds\n",
" Time in inactive batches = 1.9120E+00 seconds\n",
" Time in active batches = 1.5407E+01 seconds\n",
" Time synchronizing fission bank = 2.0000E-03 seconds\n",
" Sampling source sites = 2.0000E-03 seconds\n",
" SEND/RECV source sites = 0.0000E+00 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 1.5155E+01 seconds\n",
" Calculation Rate (inactive) = 13974.3 neutrons/second\n",
" Calculation Rate (active) = 7728.57 neutrons/second\n",
" Total time elapsed = 1.8507E+01 seconds\n",
" Calculation Rate (inactive) = 13075.3 neutrons/second\n",
" Calculation Rate (active) = 6490.56 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
@ -779,7 +779,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -890,7 +890,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -966,7 +966,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1035,7 +1035,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1111,7 +1111,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1187,7 +1187,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
"version": "2.7.10"
}
},
"nbformat": 4,

661
docs/source/pythonapi/examples/mgxs-part-ii.ipynb
View file

File diff suppressed because one or more lines are too long

227
docs/source/pythonapi/examples/mgxs-part-iii.ipynb
View file

File diff suppressed because one or more lines are too long

674
docs/source/pythonapi/examples/pandas-dataframes.ipynb
View file

File diff suppressed because one or more lines are too long

66
docs/source/pythonapi/examples/post-processing.ipynb
View file

File diff suppressed because one or more lines are too long

483
docs/source/pythonapi/examples/tally-arithmetic.ipynb
View file

@ -36,6 +36,8 @@
"import openmc\n",
"from openmc.statepoint import StatePoint\n",
"from openmc.summary import Summary\n",
"from openmc.source import Source\n",
"from openmc.stats import Box\n",
"\n",
"%matplotlib inline"
]
@ -288,7 +290,8 @@
"settings_file.particles = particles\n",
"settings_file.output = {'tallies': True, 'summary': True}\n",
"source_bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n",
"settings_file.set_source_space('box', source_bounds)\n",
"settings_file.source = Source(space=Box(\n",
" source_bounds[:3], source_bounds[3:]))\n",
"\n",
"# Export to \"settings.xml\"\n",
"settings_file.export_to_xml()"
@ -363,26 +366,7 @@
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\n",
"AAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\n",
"QYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB98LGQ4UM+6dthcAAALKSURBVGje7dpLcqQwDAbgHHE2\n",
"YeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n",
"+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\n",
"nl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n",
"/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n",
"6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\n",
"vjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\n",
"dwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\n",
"ACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\n",
"vY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n",
"+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\n",
"QBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n",
"9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n",
"8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTUtMTEtMjVUMTQ6MjA6\n",
"NTEtMDg6MDDVsKLDAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE1LTExLTI1VDE0OjIwOjUxLTA4OjAw\n",
"pO0afwAAAABJRU5ErkJggg==\n"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ABDg0ADuhPUfUAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDEtMTRUMDc6MDA6\nMTQtMDY6MDA6WZzHAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTAxLTE0VDA3OjAwOjE0LTA2OjAw\nSwQkewAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
@ -591,10 +575,9 @@
"\n",
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.0\n",
" Git SHA1: 74ffcb447521c968fb64fdaa63e40598783f2fba\n",
" Date/Time: 2015-11-25 14:20:51\n",
" MPI Processes: 1\n",
" Version: 0.7.1\n",
" Git SHA1: ea9fb637f63f9374c7436456141afa850b84acf9\n",
" Date/Time: 2016-01-14 07:00:14\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
@ -634,13 +617,13 @@
" 11/1 1.07867 1.05536 +/- 0.01277\n",
" 12/1 1.04203 1.05345 +/- 0.01096\n",
" 13/1 1.04482 1.05237 +/- 0.00955\n",
" 14/1 1.04116 1.05113 +/- 0.00852\n",
" 15/1 1.07569 1.05358 +/- 0.00800\n",
" 16/1 1.04188 1.05252 +/- 0.00732\n",
" 17/1 1.03775 1.05129 +/- 0.00679\n",
" 18/1 0.98462 1.04616 +/- 0.00808\n",
" 19/1 1.08613 1.04902 +/- 0.00801\n",
" 20/1 1.00571 1.04613 +/- 0.00800\n",
" 14/1 1.04117 1.05113 +/- 0.00852\n",
" 15/1 1.07581 1.05360 +/- 0.00801\n",
" 16/1 1.04235 1.05257 +/- 0.00731\n",
" 17/1 1.02710 1.05045 +/- 0.00701\n",
" 18/1 1.01970 1.04809 +/- 0.00687\n",
" 19/1 1.01022 1.04538 +/- 0.00691\n",
" 20/1 1.01449 1.04332 +/- 0.00675\n",
" Creating state point statepoint.20.h5...\n",
"\n",
" ===========================================================================\n",
@ -650,27 +633,27 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 7.9600E-01 seconds\n",
" Reading cross sections = 2.1200E-01 seconds\n",
" Total time in simulation = 1.8740E+01 seconds\n",
" Time in transport only = 1.8727E+01 seconds\n",
" Time in inactive batches = 2.5970E+00 seconds\n",
" Time in active batches = 1.6143E+01 seconds\n",
" Time synchronizing fission bank = 2.0000E-03 seconds\n",
" Sampling source sites = 1.0000E-03 seconds\n",
" Total time for initialization = 1.2510E+00 seconds\n",
" Reading cross sections = 9.7600E-01 seconds\n",
" Total time in simulation = 1.5844E+01 seconds\n",
" Time in transport only = 1.5834E+01 seconds\n",
" Time in inactive batches = 2.2840E+00 seconds\n",
" Time in active batches = 1.3560E+01 seconds\n",
" Time synchronizing fission bank = 3.0000E-03 seconds\n",
" Sampling source sites = 2.0000E-03 seconds\n",
" SEND/RECV source sites = 1.0000E-03 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 2.0000E-03 seconds\n",
" Total time elapsed = 1.9553E+01 seconds\n",
" Calculation Rate (inactive) = 4813.25 neutrons/second\n",
" Calculation Rate (active) = 2322.99 neutrons/second\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 1.7110E+01 seconds\n",
" Calculation Rate (inactive) = 5472.85 neutrons/second\n",
" Calculation Rate (active) = 2765.49 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.04597 +/- 0.00663\n",
" k-effective (Track-length) = 1.04613 +/- 0.00800\n",
" k-effective (Absorption) = 1.04087 +/- 0.00627\n",
" Combined k-effective = 1.04322 +/- 0.00570\n",
" k-effective (Collision) = 1.03935 +/- 0.00682\n",
" k-effective (Track-length) = 1.04332 +/- 0.00675\n",
" k-effective (Absorption) = 1.03845 +/- 0.00598\n",
" Combined k-effective = 1.04024 +/- 0.00523\n",
" Leakage Fraction = 0.00000 +/- 0.00000\n",
"\n"
]
@ -775,10 +758,10 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> total</td>\n",
" <td> (nu-fission / absorption)</td>\n",
" <td> 1.040687</td>\n",
" <td> 0.010913</td>\n",
" <td>total</td>\n",
" <td>(nu-fission / absorption)</td>\n",
" <td>1.040166</td>\n",
" <td>0.009069</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -786,7 +769,7 @@
],
"text/plain": [
" nuclide score mean std. dev.\n",
"0 total (nu-fission / absorption) 1.040687 0.010913"
"0 total (nu-fission / absorption) 1.040166 0.009069"
]
},
"execution_count": 26,
@ -836,19 +819,19 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> total</td>\n",
" <td> absorption</td>\n",
" <td> 0.959302</td>\n",
" <td> 0.010033</td>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>total</td>\n",
" <td>absorption</td>\n",
" <td>0.95938</td>\n",
" <td>0.008187</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total absorption 0.959302 0.010033"
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total absorption 0.95938 0.008187"
]
},
"execution_count": 27,
@ -896,19 +879,19 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> total</td>\n",
" <td> nu-fission</td>\n",
" <td> 1.09103</td>\n",
" <td> 0.012491</td>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>total</td>\n",
" <td>nu-fission</td>\n",
" <td>1.090899</td>\n",
" <td>0.010602</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total nu-fission 1.09103 0.012491"
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total nu-fission 1.090899 0.010602"
]
},
"execution_count": 28,
@ -958,12 +941,12 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> absorption</td>\n",
" <td> 0.803182</td>\n",
" <td> 0.008664</td>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>10000</td>\n",
" <td>total</td>\n",
" <td>absorption</td>\n",
" <td>0.803413</td>\n",
" <td>0.007031</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -971,7 +954,7 @@
],
"text/plain": [
" energy [MeV] cell nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) 10000 total absorption 0.803182 0.008664"
"0 (0.0e+00 - 6.2e-01) 10000 total absorption 0.803413 0.007031"
]
},
"execution_count": 29,
@ -1019,12 +1002,12 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> (nu-fission / absorption)</td>\n",
" <td> 1.237982</td>\n",
" <td> 0.014179</td>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>10000</td>\n",
" <td>total</td>\n",
" <td>(nu-fission / absorption)</td>\n",
" <td>1.237053</td>\n",
" <td>0.011765</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1032,10 +1015,10 @@
],
"text/plain": [
" energy [MeV] cell nuclide score mean \\\n",
"0 (0.0e+00 - 6.2e-01) 10000 total (nu-fission / absorption) 1.237982 \n",
"0 (0.0e+00 - 6.2e-01) 10000 total (nu-fission / absorption) 1.237053 \n",
"\n",
" std. dev. \n",
"0 0.014179 "
"0 0.011765 "
]
},
"execution_count": 30,
@ -1082,12 +1065,12 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> (((absorption * nu-fission) * absorption) * (n...</td>\n",
" <td> 1.040687</td>\n",
" <td> 0.022989</td>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>10000</td>\n",
" <td>total</td>\n",
" <td>(((absorption * nu-fission) * absorption) * (n...</td>\n",
" <td>1.040166</td>\n",
" <td>0.019018</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1098,7 +1081,7 @@
"0 (0.0e+00 - 6.2e-01) 10000 total \n",
"\n",
" score mean std. dev. \n",
"0 (((absorption * nu-fission) * absorption) * (n... 1.040687 0.022989 "
"0 (((absorption * nu-fission) * absorption) * (n... 1.040166 0.019018 "
]
},
"execution_count": 31,
@ -1162,75 +1145,75 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.000001</td>\n",
" <td> 8.078651e-09</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>0.000001</td>\n",
" <td>7.377419e-09</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.209990</td>\n",
" <td> 2.449396e-03</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>0.209989</td>\n",
" <td>2.303838e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.356117</td>\n",
" <td> 4.364366e-03</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>0.356420</td>\n",
" <td>3.951669e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.005555</td>\n",
" <td> 6.495710e-05</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>0.005555</td>\n",
" <td>6.101004e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.007190</td>\n",
" <td> 7.596666e-05</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>0.007155</td>\n",
" <td>8.053460e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.227843</td>\n",
" <td> 1.024510e-03</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>0.227770</td>\n",
" <td>1.079289e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.008086</td>\n",
" <td> 6.251590e-05</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>0.008067</td>\n",
" <td>5.254797e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.003365</td>\n",
" <td> 1.646663e-05</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>0.003367</td>\n",
" <td>1.647058e-05</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1239,23 +1222,23 @@
"text/plain": [
" cell energy [MeV] nuclide score mean \\\n",
"0 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (nu-fission / flux) 0.000001 \n",
"1 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (scatter / flux) 0.209990 \n",
"2 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (nu-fission / flux) 0.356117 \n",
"1 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (scatter / flux) 0.209989 \n",
"2 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (nu-fission / flux) 0.356420 \n",
"3 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (scatter / flux) 0.005555 \n",
"4 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (nu-fission / flux) 0.007190 \n",
"5 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (scatter / flux) 0.227843 \n",
"6 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (nu-fission / flux) 0.008086 \n",
"7 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (scatter / flux) 0.003365 \n",
"4 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (nu-fission / flux) 0.007155 \n",
"5 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (scatter / flux) 0.227770 \n",
"6 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (nu-fission / flux) 0.008067 \n",
"7 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (scatter / flux) 0.003367 \n",
"\n",
" std. dev. \n",
"0 8.078651e-09 \n",
"1 2.449396e-03 \n",
"2 4.364366e-03 \n",
"3 6.495710e-05 \n",
"4 7.596666e-05 \n",
"5 1.024510e-03 \n",
"6 6.251590e-05 \n",
"7 1.646663e-05 "
"0 7.377419e-09 \n",
"1 2.303838e-03 \n",
"2 3.951669e-03 \n",
"3 6.101004e-05 \n",
"4 8.053460e-05 \n",
"5 1.079289e-03 \n",
"6 5.254797e-05 \n",
"7 1.647058e-05 "
]
},
"execution_count": 33,
@ -1286,11 +1269,11 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 6.65302296e-07]\n",
" [ 3.56116716e-01]]\n",
"[[[ 6.65702880e-07]\n",
" [ 3.56420449e-01]]\n",
"\n",
" [[ 7.19004460e-03]\n",
" [ 8.08598751e-03]]]\n"
" [[ 7.15488656e-03]\n",
" [ 8.06673774e-03]]]\n"
]
}
],
@ -1318,9 +1301,9 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.00555516]]\n",
"[[[ 0.00555533]]\n",
"\n",
" [[ 0.00336498]]]\n"
" [[ 0.0033668 ]]]\n"
]
}
],
@ -1342,8 +1325,8 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.22784316]\n",
" [ 0.00336498]]]\n"
"[[[ 0.22777006]\n",
" [ 0.0033668 ]]]\n"
]
}
],
@ -1388,39 +1371,39 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> U-238</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.000002</td>\n",
" <td> 1.450189e-08</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>U-238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.000002</td>\n",
" <td>1.283958e-08</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> U-235</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.870882</td>\n",
" <td> 7.895515e-03</td>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>U-235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.868553</td>\n",
" <td>6.880390e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> U-238</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.082484</td>\n",
" <td> 8.253437e-04</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>U-238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.082149</td>\n",
" <td>8.837250e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> U-235</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.092762</td>\n",
" <td> 6.444580e-04</td>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>U-235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.092618</td>\n",
" <td>5.195308e-04</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1428,10 +1411,10 @@
],
"text/plain": [
" cell energy [MeV] nuclide score mean std. dev.\n",
"0 10000 (0.0e+00 - 6.3e-07) U-238 nu-fission 0.000002 1.450189e-08\n",
"1 10000 (0.0e+00 - 6.3e-07) U-235 nu-fission 0.870882 7.895515e-03\n",
"2 10000 (6.3e-07 - 2.0e+01) U-238 nu-fission 0.082484 8.253437e-04\n",
"3 10000 (6.3e-07 - 2.0e+01) U-235 nu-fission 0.092762 6.444580e-04"
"0 10000 (0.0e+00 - 6.3e-07) U-238 nu-fission 0.000002 1.283958e-08\n",
"1 10000 (0.0e+00 - 6.3e-07) U-235 nu-fission 0.868553 6.880390e-03\n",
"2 10000 (6.3e-07 - 2.0e+01) U-238 nu-fission 0.082149 8.837250e-04\n",
"3 10000 (6.3e-07 - 2.0e+01) U-235 nu-fission 0.092618 5.195308e-04"
]
},
"execution_count": 37,
@ -1471,84 +1454,84 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 10002</td>\n",
" <td> (1.0e-08 - 1.1e-07)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 4.630154</td>\n",
" <td> 0.044512</td>\n",
" <td>10002</td>\n",
" <td>(1.0e-08 - 1.1e-07)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>4.619398</td>\n",
" <td>0.040124</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 10002</td>\n",
" <td> (1.1e-07 - 1.2e-06)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.042984</td>\n",
" <td> 0.011429</td>\n",
" <td>10002</td>\n",
" <td>(1.1e-07 - 1.2e-06)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.030757</td>\n",
" <td>0.011239</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td> 10002</td>\n",
" <td> (1.2e-06 - 1.3e-05)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 1.657517</td>\n",
" <td> 0.008617</td>\n",
" <td>10002</td>\n",
" <td>(1.2e-06 - 1.3e-05)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>1.658488</td>\n",
" <td>0.009777</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td> 10002</td>\n",
" <td> (1.3e-05 - 1.4e-04)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 1.863326</td>\n",
" <td> 0.008848</td>\n",
" <td>10002</td>\n",
" <td>(1.3e-05 - 1.4e-04)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>1.853002</td>\n",
" <td>0.007378</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td> 10002</td>\n",
" <td> (1.4e-04 - 1.5e-03)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.043916</td>\n",
" <td> 0.014195</td>\n",
" <td>10002</td>\n",
" <td>(1.4e-04 - 1.5e-03)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.050773</td>\n",
" <td>0.012484</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td> 10002</td>\n",
" <td> (1.5e-03 - 1.6e-02)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.134458</td>\n",
" <td> 0.007561</td>\n",
" <td>10002</td>\n",
" <td>(1.5e-03 - 1.6e-02)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.131759</td>\n",
" <td>0.007821</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td> 10002</td>\n",
" <td> (1.6e-02 - 1.7e-01)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.209947</td>\n",
" <td> 0.013848</td>\n",
" <td>10002</td>\n",
" <td>(1.6e-02 - 1.7e-01)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.213710</td>\n",
" <td>0.015159</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td> 10002</td>\n",
" <td> (1.7e-01 - 1.9e+00)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.006967</td>\n",
" <td> 0.009368</td>\n",
" <td>10002</td>\n",
" <td>(1.7e-01 - 1.9e+00)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.011925</td>\n",
" <td>0.009406</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td> 10002</td>\n",
" <td> (1.9e+00 - 2.0e+01)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 0.373895</td>\n",
" <td> 0.002964</td>\n",
" <td>10002</td>\n",
" <td>(1.9e+00 - 2.0e+01)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>0.371280</td>\n",
" <td>0.003949</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1556,15 +1539,15 @@
],
"text/plain": [
" cell energy [MeV] nuclide score mean std. dev.\n",
"0 10002 (1.0e-08 - 1.1e-07) H-1 scatter 4.630154 0.044512\n",
"1 10002 (1.1e-07 - 1.2e-06) H-1 scatter 2.042984 0.011429\n",
"2 10002 (1.2e-06 - 1.3e-05) H-1 scatter 1.657517 0.008617\n",
"3 10002 (1.3e-05 - 1.4e-04) H-1 scatter 1.863326 0.008848\n",
"4 10002 (1.4e-04 - 1.5e-03) H-1 scatter 2.043916 0.014195\n",
"5 10002 (1.5e-03 - 1.6e-02) H-1 scatter 2.134458 0.007561\n",
"6 10002 (1.6e-02 - 1.7e-01) H-1 scatter 2.209947 0.013848\n",
"7 10002 (1.7e-01 - 1.9e+00) H-1 scatter 2.006967 0.009368\n",
"8 10002 (1.9e+00 - 2.0e+01) H-1 scatter 0.373895 0.002964"
"0 10002 (1.0e-08 - 1.1e-07) H-1 scatter 4.619398 0.040124\n",
"1 10002 (1.1e-07 - 1.2e-06) H-1 scatter 2.030757 0.011239\n",
"2 10002 (1.2e-06 - 1.3e-05) H-1 scatter 1.658488 0.009777\n",
"3 10002 (1.3e-05 - 1.4e-04) H-1 scatter 1.853002 0.007378\n",
"4 10002 (1.4e-04 - 1.5e-03) H-1 scatter 2.050773 0.012484\n",
"5 10002 (1.5e-03 - 1.6e-02) H-1 scatter 2.131759 0.007821\n",
"6 10002 (1.6e-02 - 1.7e-01) H-1 scatter 2.213710 0.015159\n",
"7 10002 (1.7e-01 - 1.9e+00) H-1 scatter 2.011925 0.009406\n",
"8 10002 (1.9e+00 - 2.0e+01) H-1 scatter 0.371280 0.003949"
]
},
"execution_count": 38,

2
docs/source/pythonapi/index.rst
View file

@ -35,6 +35,8 @@ on a given module or class.
opencg_compatible
plots
settings
source
stats
surface
tallies
trigger

8
docs/source/pythonapi/source.rst Normal file
View file

@ -0,0 +1,8 @@
.. _pythonapi_source:
======
Source
======
.. automodule:: openmc.source
:members:

58
docs/source/pythonapi/stats.rst Normal file
View file

@ -0,0 +1,58 @@
.. _pythonapi_stats:
=====================
Statistical Functions
=====================
----------------------------
Summary of Available Classes
----------------------------
Univariate Probability Distributions
------------------------------------
.. currentmodule:: openmc.stats.univariate
.. autosummary::
Univariate
Discrete
Uniform
Maxwell
Watt
Tabular
Angular Distributions
---------------------
.. currentmodule:: openmc.stats.multivariate
.. autosummary::
UnitSphere
PolarAzimuthal
Isotropic
Monodirectional
Spatial Distributions
---------------------
.. autosummary::
Spatial
CartesianIndependent
Box
Point
Univariate Probability Distributions
------------------------------------
.. automodule:: openmc.stats.univariate
:members:
Multivariate Probability Distributions
--------------------------------------
.. automodule:: openmc.stats.multivariate
:members:

166
docs/source/usersguide/input.rst
View file

@ -424,9 +424,17 @@ pseudo-random number generator.
The ``source`` element gives information on an external source distribution to
be used either as the source for a fixed source calculation or the initial
source guess for criticality calculations. It takes the following
source guess for criticality calculations. Multiple ``<source>`` elements may be
specified to define different source distributions. Each one takes the following
attributes/sub-elements:
:strength:
The strength of the source. If multiple sources are present, the source
strength indicates the relative probability of choosing one source over the
other.
*Default*: 1.0
:file:
If this attribute is given, it indicates that the source is to be read from
a binary source file whose path is given by the value of this element. Note,
@ -440,13 +448,13 @@ attributes/sub-elements:
has the following attributes:
:type:
The type of spatial distribution. Valid options are "box", "fission", and
"point". A "box" spatial distribution has coordinates sampled uniformly in
a parallelepiped. A "fission" spatial distribution samples locations from
a "box" distribution but only locations in fissionable materials are
accepted. A "point" spatial distribution has coordinates specified by a
triplet.
The type of spatial distribution. Valid options are "box", "fission",
"point", and "cartesian". A "box" spatial distribution has coordinates
sampled uniformly in a parallelepiped. A "fission" spatial distribution
samples locations from a "box" distribution but only locations in
fissionable materials are accepted. A "point" spatial distribution has
coordinates specified by a triplet. An "cartesian" spatial distribution
specifies independent distributions of x-, y-, and z-coordinates.
*Default*: None
@ -459,67 +467,123 @@ attributes/sub-elements:
For a "point" spatial distribution, ``parameters`` should be given as
three real numbers which specify the (x,y,z) location of an isotropic
point source
point source.
For an "cartesian" distribution, no parameters are specified. Instead,
the ``x``, ``y``, and ``z`` elements must be specified.
*Default*: None
:x:
For an "cartesian" distribution, this element specifies the distribution
of x-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
:ref:`univariate`).
:y:
For an "cartesian" distribution, this element specifies the distribution
of y-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
:ref:`univariate`).
:z:
For an "cartesian" distribution, this element specifies the distribution
of z-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
:ref:`univariate`).
:angle:
An element specifying the angular distribution of source sites. This element
has the following attributes:
:type:
The type of angular distribution. Valid options are "isotropic" and
"monodirectional". The angle of the particle emitted from a source site is
isotropic if the "isotropic" option is given. The angle of the particle
emitted from a source site is the direction specified in the <parameters>
attribute if "monodirectional" option is given.
The type of angular distribution. Valid options are "isotropic",
"monodirectional", and "mu-phi". The angle of the particle emitted from a
source site is isotropic if the "isotropic" option is given. The angle of
the particle emitted from a source site is the direction specified in the
``reference_uvw`` element/attribute if "monodirectional" option is
given. The "mu-phi" option produces directions with the cosine of the
polar angle and the azimuthal angle explicitly specified.
*Default*: isotropic
:parameters:
For an "isotropic" angular distribution, ``parameters`` should not be
specified.
:reference_uvw:
The direction from which the polar angle is measured. Represented by the
x-, y-, and z-components of a unit vector. For a monodirectional
distribution, this defines the direction of all sampled particles.
For a "monodirectional" angular distribution, ``parameters`` should be
given as three real numbers which specify the angular cosines with respect
to each axis.
:mu:
An element specifying the distribution of the cosine of the polar
angle. Only relevant when the type is "mu-phi". The necessary
sub-elements/attributes are those of a univariate probability distribution
(see the description in :ref:`univariate`).
*Default*: None
:phi:
An element specifying the distribution of the azimuthal angle. Only
relevant when the type is "mu-phi". The necessary sub-elements/attributes
are those of a univariate probability distribution (see the description in
:ref:`univariate`).
:energy:
An element specifying the energy distribution of source sites. This element
has the following attributes:
:type:
The type of energy distribution. Valid options are "monoenergetic",
"watt", and "maxwell". The "monoenergetic" option produces source sites at
a single energy. The "watt" option produces source sites whose energy is
sampled from a Watt fission spectrum. The "maxwell" option produce source
sites whose energy is sampled from a Maxwell fission spectrum.
*Default*: watt
:parameters:
For a "monoenergetic" energy distribution, ``parameters`` should be
given as the energy in MeV of the source sites.
For a "watt" energy distribution, ``parameters`` should be given as two
real numbers :math:`a` and :math:`b` that parameterize the distribution
:math:`p(E) dE = c e^{-E/a} \sinh \sqrt{b \, E} dE`.
For a "maxwell" energy distribution, ``parameters`` should be given as one
real number :math:`a` that parameterizes the distribution :math:`p(E) dE =
c E e^{-E/a} dE`.
*Default*: 0.988 2.249
An element specifying the energy distribution of source sites. The necessary
sub-elements/attributes are those of a univariate probability distribution
(see the description in :ref:`univariate`).
:write_initial:
An element specifying whether to write out the initial source bank used at
the beginning of the first batch. The output file is named
"initial_source.binary(h5)"
"initial_source.h5"
*Default*: false
*Default*: false
.. _univariate:
Univariate Probability Distributions
++++++++++++++++++++++++++++++++++++
Various components of a source distribution involve probability distributions of
a single random variable, e.g. the distribution of the energy, the distribution
of the polar angle, and the distribution of x-coordinates. Each of these
components supports the same syntax with an element whose tag signifies the
variable and whose sub-elements/attributes are as follows:
:type:
The type of the distribution. Valid options are "uniform", "discrete",
"tabular", "maxwell", and "watt". The "uniform" option produces variates
sampled from a uniform distribution over a finite interval. The "discrete"
option produces random variates that can assume a finite number of values
(i.e., a distribution characterized by a probability mass function). The
"tabular" option produces random variates sampled from a tabulated
distribution where the density function is either a histogram or
linearly-interpolated between tabulated points. The "watt" option produces
random variates is sampled from a Watt fission spectrum (only used for
energies). The "maxwell" option produce variates sampled from a Maxwell
fission spectrum (only used for energies).
*Default*: None
:parameters:
For a "uniform" distribution, ``parameters`` should be given as two real
numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over
which random variates are sampled.
For a "discrete" or "tabular" distribution, ``parameters`` provides the
:math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
points are given first followed by corresponding :math:`p` points.
For a "watt" distribution, ``parameters`` should be given as two real numbers
:math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
e^{-x/a} \sinh \sqrt{b \, x} dx`.
For a "maxwell" distribution, ``parameters`` should be given as one real
number :math:`a` that parameterizes the distribution :math:`p(x) dx = c x
e^{-x/a} dx`.
:interpolation:
For a "tabular" distribution, ``interpolation`` can be set to "histogram" or
"linear-linear" thereby specifying how tabular points are to be interpolated.
*Default*: histogram
``<state_point>`` Element
-------------------------
@ -891,7 +955,9 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
:material:
The ``id`` of the material that this cell contains. If the cell should
contain no material, this can also be set to "void".
contain no material, this can also be set to "void". A list of materials
can be specified for the "distributed material" feature. This will give each
unique instance of the cell its own material.
.. note:: If a material is specified, no fill should be given.

4
docs/source/usersguide/output/statepoint.rst
View file

@ -185,10 +185,6 @@ if run_mode == 'k-eigenvalue':
Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn',
'surface', 'mesh', 'energy', 'energyout', or 'distribcell'.
**/tallies/tally <uid>/filter <j>/offset** (*int*)
Filter offset (used for distribcell filter).
**/tallies/tally <uid>/filter <j>/n_bins** (*int*)
Number of bins for the j-th filter.

17
docs/source/usersguide/output/summary.rst
View file

@ -91,10 +91,12 @@ The current revision of the summary file format is 1.
Type of fill for the cell. Can be 'normal', 'universe', or 'lattice'.
**/geometry/cells/cell <uid>/material** (*int*)
**/geometry/cells/cell <uid>/material** (*int* or *int[]*)
Unique ID of the material assigned to the cell. This dataset is present only
if fill_type is set to 'normal'.
Unique ID of the material(s) assigned to the cell. This dataset is present
only if fill_type is set to 'normal'. The value '-1' signifies void
material. The data is an array if the cell uses distributed materials,
otherwise it is a scalar.
**/geometry/cells/cell <uid>/offset** (*int[]*)
@ -121,6 +123,10 @@ The current revision of the summary file format is 1.
Region specification for the cell.
**/geometry/cells/cell <uid>/distribcell_index** (*int*)
Index of this cell in distribcell filter arrays.
**/geometry/surfaces/surface <uid>/index** (*int*)
Index in surfaces array used internally in OpenMC.
@ -306,6 +312,11 @@ The current revision of the summary file format is 1.
than the number of user-specified scores since each score might have
multiple scoring bins, e.g., scatter-PN.
**/tallies/tally <uid>/moment_orders** (*char[][]*)
Tallying moment orders for Legendre and spherical harmonic tally expansions
(*e.g.*, 'P2', 'Y1,2', etc.).
**/tallies/tally <uid>/score_bins** (*char[][]*)
Scoring bins for the tally.

5
docs/sphinxext/notebook_sphinxext.py
View file

@ -12,10 +12,7 @@ from docutils.parsers.rst import directives, roles, states
from docutils.parsers.rst.roles import set_classes
from docutils.transforms import misc
try:
from IPython.nbconver.exporters import html
except ImportError:
from IPython.nbconvert import html
from nbconvert import html
class Notebook(Directive):

4
examples/python/basic/build-xml.py
View file

@ -1,4 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -92,7 +94,7 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-4, -4, -4, 4, 4, 4])
settings_file.source = Source(space=Box([-4, -4, -4], [4, 4, 4]))
settings_file.export_to_xml()

4
examples/python/boxes/build-xml.py
View file

@ -1,6 +1,8 @@
import numpy as np
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -117,7 +119,7 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', np.concatenate(outer_cube.bounding_box))
settings_file.source = Source(space=Box(*outer_cube.bounding_box))
settings_file.export_to_xml()
###############################################################################

6
examples/python/lattice/hexagonal/build-xml.py
View file

@ -1,5 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -125,7 +126,8 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
settings_file.keff_trigger = {'type' : 'std_dev', 'threshold' : 5E-4}
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100

6
examples/python/lattice/nested/build-xml.py
View file

@ -1,5 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -136,7 +137,8 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
settings_file.export_to_xml()

5
examples/python/lattice/simple/build-xml.py
View file

@ -1,4 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -125,7 +127,8 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100
settings_file.export_to_xml()

6
examples/python/pincell/build-xml.py
View file

@ -1,4 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -168,8 +170,8 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-0.62992, -0.62992, -1, \
0.62992, 0.62992, 1])
settings_file.source = Source(space=Box(
[-0.62992, -0.62992, -1], [0.62992, 0.62992, 1]))
settings_file.entropy_lower_left = [-0.39218, -0.39218, -1.e50]
settings_file.entropy_upper_right = [0.39218, 0.39218, 1.e50]
settings_file.entropy_dimension = [10, 10, 1]

4
examples/python/reflective/build-xml.py
View file

@ -1,6 +1,8 @@
import numpy as np
import openmc
from openmc.stats import Box
from openmc.source import Source
###############################################################################
# Simulation Input File Parameters
@ -84,5 +86,5 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', np.concatenate(cell.region.bounding_box))
settings_file.source = Source(space=Box(*cell.region.bounding_box))
settings_file.export_to_xml()

409
openmc/aggregate.py Normal file
View file

@ -0,0 +1,409 @@
import sys
from numbers import Integral
import numpy as np
from openmc import Filter, Nuclide
from openmc.cross import CrossScore, CrossNuclide, CrossFilter
from openmc.filter import _FILTER_TYPES
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
# Acceptable tally aggregation operations
_TALLY_AGGREGATE_OPS = ['sum', 'mean']
class AggregateScore(object):
"""A special-purpose tally score used to encapsulate an aggregate of a
subset or all of tally's scores for tally aggregation.
Parameters
----------
scores : Iterable of str or CrossScore
The scores included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally's scores with this AggregateScore
Attributes
----------
scores : Iterable of str or CrossScore
The scores included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally's scores with this AggregateScore
"""
def __init__(self, scores=None, aggregate_op=None):
self._scores = None
self._aggregate_op = None
if scores is not None:
self.scores = scores
if aggregate_op is not None:
self.aggregate_op = aggregate_op
def __hash__(self):
return hash(repr(self))
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._scores = self.scores
clone._aggregate_op = self.aggregate_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
string = ', '.join(map(str, self.scores))
string = '{0}({1})'.format(self.aggregate_op, string)
return string
@property
def scores(self):
return self._scores
@property
def aggregate_op(self):
return self._aggregate_op
@scores.setter
def scores(self, scores):
cv.check_iterable_type('scores', scores, basestring)
self._scores = scores
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, (basestring, CrossScore))
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op
class AggregateNuclide(object):
"""A special-purpose tally nuclide used to encapsulate an aggregate of a
subset or all of tally's nuclides for tally aggregation.
Parameters
----------
nuclides : Iterable of str or Nuclide or CrossNuclide
The nuclides included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally's nuclides with this AggregateNuclide
Attributes
----------
nuclides : Iterable of str or Nuclide or CrossNuclide
The nuclides included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally's nuclides with this AggregateNuclide
"""
def __init__(self, nuclides=None, aggregate_op=None):
self._nuclides = None
self._aggregate_op = None
if nuclides is not None:
self.nuclides = nuclides
if aggregate_op is not None:
self.aggregate_op = aggregate_op
def __hash__(self):
return hash(repr(self))
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._nuclides = self.nuclides
clone._aggregate_op = self._aggregate_op
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
def __repr__(self):
# Append each nuclide in the aggregate to the string
string = '{0}('.format(self.aggregate_op)
names = [nuclide.name if isinstance(nuclide, Nuclide) else str(nuclide)
for nuclide in self.nuclides]
string += ', '.join(map(str, names)) + ')'
return string
@property
def nuclides(self):
return self._nuclides
@property
def aggregate_op(self):
return self._aggregate_op
@nuclides.setter
def nuclides(self, nuclides):
cv.check_iterable_type('nuclides', nuclides,
(basestring, Nuclide, CrossNuclide))
self._nuclides = nuclides
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, basestring)
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op
class AggregateFilter(object):
"""A special-purpose tally filter used to encapsulate an aggregate of a
subset or all of a tally filter's bins for tally aggregation.
Parameters
----------
aggregate_filter : Filter or CrossFilter
The filter included in the aggregation
bins : Iterable of tuple
The filter bins included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally filter's bins with this AggregateFilter
Attributes
----------
type : str
The type of the aggregatefilter (e.g., 'sum(energy)', 'sum(cell)')
aggregate_filter : filter
The filter included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
to aggregate across a tally filter's bins with this AggregateFilter
bins : Iterable of tuple
The filter bins included in the aggregation
num_bins : Integral
The number of filter bins (always 1 if aggregate_filter is defined)
stride : Integral
The number of filter, nuclide and score bins within each of this
aggregatefilter's bins.
"""
def __init__(self, aggregate_filter=None, bins=None, aggregate_op=None):
self._type = '{0}({1})'.format(aggregate_op, aggregate_filter.type)
self._bins = None
self._stride = None
self._aggregate_filter = None
self._aggregate_op = None
if aggregate_filter is not None:
self.aggregate_filter = aggregate_filter
if bins is not None:
self.bins = bins
if aggregate_op is not None:
self.aggregate_op = aggregate_op
def __hash__(self):
return hash((self.type, self.bins, self.aggregate_op))
def __eq__(self, other):
return str(other) == str(self)
def __ne__(self, other):
return not self == other
def __repr__(self):
string = 'AggregateFilter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
return string
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is the first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._type = self.type
clone._aggregate_filter = self.aggregate_filter
clone._aggregate_op = self.aggregate_op
clone._bins = self._bins
clone._stride = self.stride
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def aggregate_filter(self):
return self._aggregate_filter
@property
def aggregate_op(self):
return self._aggregate_op
@property
def type(self):
return self._type
@property
def bins(self):
return self._bins
@property
def num_bins(self):
return 1 if self.aggregate_filter else 0
@property
def stride(self):
return self._stride
@type.setter
def type(self, filter_type):
if filter_type not in _FILTER_TYPES.values():
msg = 'Unable to set AggregateFilter type to "{0}" since it ' \
'is not one of the supported types'.format(filter_type)
raise ValueError(msg)
self._type = filter_type
@aggregate_filter.setter
def aggregate_filter(self, aggregate_filter):
cv.check_type('aggregate_filter', aggregate_filter, (Filter, CrossFilter))
self._aggregate_filter = aggregate_filter
@bins.setter
def bins(self, bins):
cv.check_iterable_type('bins', bins, (Integral, tuple))
self._bins = bins
@aggregate_op.setter
def aggregate_op(self, aggregate_op):
cv.check_type('aggregate_op', aggregate_op, basestring)
cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
self._aggregate_op = aggregate_op
@stride.setter
def stride(self, stride):
self._stride = stride
def get_bin_index(self, filter_bin):
"""Returns the index in the AggregateFilter for some bin.
Parameters
----------
filter_bin : Integral or tuple of Real
A tuple of value(s) corresponding to the bin of interest in
the aggregated filter. The bin is the integer ID for 'material',
'surface', 'cell', 'cellborn', and 'universe' Filters. The bin
is the integer cell instance ID for 'distribcell' Filters. The
bin is a 2-tuple of floats for 'energy' and 'energyout' filters
corresponding to the energy boundaries of the bin of interest.
The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to
the mesh cell of interest.
Returns
-------
filter_index : Integral
The index in the Tally data array for this filter bin. For an
AggregateTally the filter bin index is always unity.
Raises
------
ValueError
When the filter_bin is not part of the aggregated filter's bins
"""
if filter_bin not in self.bins:
msg = 'Unable to get the bin index for AggregateFilter since ' \
'"{0}" is not one of the bins'.format(filter_bin)
raise ValueError(msg)
else:
return 0
def get_pandas_dataframe(self, datasize, summary=None):
"""Builds a Pandas DataFrame for the AggregateFilter's bins.
This method constructs a Pandas DataFrame object for the AggregateFilter
with columns annotated by filter bin information. This is a helper
method for the Tally.get_pandas_dataframe(...) method.
Parameters
----------
datasize : Integral
The total number of bins in the tally corresponding to this filter
summary : None or Summary
An optional Summary object to be used to construct columns for
distribcell tally filters (default is None). NOTE: This parameter
is not used by the AggregateFilter and simply mirrors the method
signature for the CrossFilter.
Returns
-------
pandas.DataFrame
A Pandas DataFrame with columns of strings that characterize the
aggregatefilter's bins. Each entry in the DataFrame will include
one or more aggregation operations used to construct the
aggregatefilter's bins. The number of rows in the DataFrame is the
same as the total number of bins in the corresponding tally, with
the filter bins appropriately tiled to map to the corresponding
tally bins.
See also
--------
Tally.get_pandas_dataframe(), Filter.get_pandas_dataframe(),
CrossFilter.get_pandas_dataframe()
"""
import pandas as pd
# Construct a sring representing the filter aggregation
aggregate_bin = '{0}('.format(self.aggregate_op)
aggregate_bin += ', '.join(map(str, self.bins)) + ')'
# Construct NumPy array of bin repeated for each element in dataframe
aggregate_bin_array = np.array([aggregate_bin])
aggregate_bin_array = np.repeat(aggregate_bin_array, datasize)
# Construct Pandas DataFrame for the AggregateFilter
df = pd.DataFrame({self.aggregate_filter.type: aggregate_bin_array})
return df

17
openmc/cross.py
View file

@ -34,7 +34,7 @@ class CrossScore(object):
The right score in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's scores with this CrossNuclide
combine two tally's scores with this CrossScore
"""
@ -245,6 +245,8 @@ class CrossFilter(object):
Attributes
----------
type : str
The type of the crossfilter (e.g., 'energy / energy')
left_filter : Filter or CrossFilter
The left filter in the outer product
right_filter : Filter or CrossFilter
@ -252,6 +254,14 @@ class CrossFilter(object):
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
combine two tally's filter bins with this CrossFilter
bins : dict of Iterable
A dictionary of the bins from each filter keyed by the types of the
left / right filters
num_bins : Integral
The number of filter bins (always 1 if aggregate_filter is defined)
stride : Integral
The number of filter, nuclide and score bins within each of this
crossfilter's bins.
"""
@ -310,7 +320,6 @@ class CrossFilter(object):
clone._binary_op = self.binary_op
clone._type = self.type
clone._bins = self._bins
clone._num_bins = self.num_bins
clone._stride = self.stride
memo[id(self)] = clone
@ -355,8 +364,8 @@ class CrossFilter(object):
@type.setter
def type(self, filter_type):
if filter_type not in _FILTER_TYPES.values():
msg = 'Unable to set Filter type to "{0}" since it is not one ' \
'of the supported types'.format(filter_type)
msg = 'Unable to set CrossFilter type to "{0}" since it ' \
'is not one of the supported types'.format(filter_type)
raise ValueError(msg)
self._type = filter_type

38
openmc/filter.py
View file

@ -35,15 +35,13 @@ class Filter(object):
Attributes
----------
type : str
The type of the tally filter.
bins : Integral or Iterable of Integral or Iterable of Real
The type of the tally filter
bins : Integral or Iterable of Real
The bins for the filter
num_bins : Integral
The number of filter bins
mesh : Mesh or None
A Mesh object for 'mesh' type filters.
offset : Integral
A value used to index tally bins for 'distribcell' tallies.
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
@ -57,7 +55,6 @@ class Filter(object):
self._num_bins = 0
self._bins = None
self._mesh = None
self._offset = -1
self._stride = None
if type is not None:
@ -93,7 +90,6 @@ class Filter(object):
clone._bins = copy.deepcopy(self.bins, memo)
clone._num_bins = self.num_bins
clone._mesh = copy.deepcopy(self.mesh, memo)
clone._offset = self.offset
clone._stride = self.stride
memo[id(self)] = clone
@ -108,7 +104,6 @@ class Filter(object):
string = 'Filter\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self.offset)
return string
@property
@ -134,10 +129,6 @@ class Filter(object):
def mesh(self):
return self._mesh
@property
def offset(self):
return self._offset
@property
def stride(self):
return self._stride
@ -226,11 +217,6 @@ class Filter(object):
self.type = 'mesh'
self.bins = self.mesh.id
@offset.setter
def offset(self, offset):
cv.check_type('filter offset', offset, Integral)
self._offset = offset
@stride.setter
def stride(self, stride):
cv.check_type('filter stride', stride, Integral)
@ -241,12 +227,12 @@ class Filter(object):
self._stride = stride
def can_merge(self, filter):
def can_merge(self, other):
"""Determine if filter can be merged with another.
Parameters
----------
filter : Filter
other : Filter
Filter to compare with
Returns
@ -256,11 +242,11 @@ class Filter(object):
"""
if not isinstance(filter, Filter):
if not isinstance(other, Filter):
return False
# Filters must be of the same type
elif self.type != filter.type:
elif self.type != other.type:
return False
# Distribcell filters cannot have more than one bin
@ -278,12 +264,12 @@ class Filter(object):
else:
return True
def merge(self, filter):
def merge(self, other):
"""Merge this filter with another.
Parameters
----------
filter : Filter
other : Filter
Filter to merge with
Returns
@ -293,16 +279,16 @@ class Filter(object):
"""
if not self.can_merge(filter):
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(self.type, filter.type)
'filters'.format(self.type, other.type)
raise ValueError(msg)
# Create deep copy of filter to return as merged filter
merged_filter = copy.deepcopy(self)
# Merge unique filter bins
merged_bins = list(set(np.concatenate((self.bins, filter.bins))))
merged_bins = list(set(np.concatenate((self.bins, other.bins))))
merged_filter.bins = merged_bins
merged_filter.num_bins = len(merged_bins)
@ -623,7 +609,7 @@ class Filter(object):
# If this region is in Cell corresponding to the
# distribcell filter bin, store it in dictionary
if cell_id == self.bins[0]:
offset = openmc_geometry.get_offset(path, self.offset)
offset = openmc_geometry.get_cell_instance(path)
offsets_to_coords[offset] = coords
# Each distribcell offset is a DataFrame bin

38
openmc/geometry.py
View file

@ -1,4 +1,4 @@
from collections import OrderedDict
from collections import Iterable, OrderedDict
from xml.etree import ElementTree as ET
import openmc
@ -42,10 +42,10 @@ class Geometry(object):
self._root_universe = root_universe
def get_offset(self, path, filter_offset):
"""Returns the corresponding location in the results array for a given path and
filter number. This is primarily intended to post-processing result when
a distribcell filter is used.
def get_cell_instance(self, path):
"""Return the instance number for the final cell in a geometry path.
The instance is an index into tally distribcell filter arrays.
Parameters
----------
@ -55,24 +55,31 @@ class Geometry(object):
lattice passed through. For the case of the lattice, a tuple should
be provided to indicate which coordinates in the lattice should be
entered. This should be in the form: (lat_id, i_x, i_y, i_z)
filter_offset : int
An integer that specifies which offset map the filter is using
Returns
-------
offset : int
Location in the results array for the path and filter
instance : int
Index in tally results array for distribcell filters
"""
# Find the distribcell index of the cell.
cells = self.get_all_cells()
if path[-1] in cells:
distribcell_index = cells[path[-1]].distribcell_index
else:
raise RuntimeError('Could not find cell {} specified in a \
distribcell filter'.format(path[-1]))
# Return memoize'd offset if possible
if (path, filter_offset) in self._offsets:
offset = self._offsets[(path, filter_offset)]
if (path, distribcell_index) in self._offsets:
offset = self._offsets[(path, distribcell_index)]
# Begin recursive call to compute offset starting with the base Universe
else:
offset = self._root_universe.get_offset(path, filter_offset)
self._offsets[(path, filter_offset)] = offset
offset = self._root_universe.get_cell_instance(path,
distribcell_index)
self._offsets[(path, distribcell_index)] = offset
# Return the final offset
return offset
@ -133,7 +140,10 @@ class Geometry(object):
materials = set()
for cell in material_cells:
materials.add(cell._fill)
if isinstance(cell.fill, Iterable):
for m in cell.fill: materials.add(m)
else:
materials.add(cell.fill)
materials = list(materials)
materials.sort(key=lambda x: x.id)

32
openmc/mgxs/library.py
View file

@ -73,6 +73,9 @@ class Library(object):
name : str, optional
Name of the multi-group cross section library. Used as a label to
identify tallies in OpenMC 'tallies.xml' file.
sparse : bool
Whether or not the Library's tallies use SciPy's LIL sparse matrix
format for compressed data storage
"""
@ -92,6 +95,7 @@ class Library(object):
self._all_mgxs = OrderedDict()
self._sp_filename = None
self._keff = None
self._sparse = False
self.name = name
self.openmc_geometry = openmc_geometry
@ -119,6 +123,7 @@ class Library(object):
clone._all_mgxs = self.all_mgxs
clone._sp_filename = self._sp_filename
clone._keff = self._keff
clone._sparse = self.sparse
clone._all_mgxs = OrderedDict()
for domain in self.domains:
@ -208,6 +213,10 @@ class Library(object):
def keff(self):
return self._keff
@property
def sparse(self):
return self._sparse
@openmc_geometry.setter
def openmc_geometry(self, openmc_geometry):
cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
@ -286,6 +295,28 @@ class Library(object):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
self._tally_trigger = tally_trigger
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
# Sparsify or densify each MGXS in the Library
for domain in self.domains:
for mgxs_type in self.mgxs_types:
mgxs = self.get_mgxs(domain, mgxs_type)
mgxs.sparse = self.sparse
self._sparse = sparse
def build_library(self):
"""Initialize MGXS objects in each domain and for each reaction type
in the library.
@ -381,6 +412,7 @@ class Library(object):
for mgxs_type in self.mgxs_types:
mgxs = self.get_mgxs(domain, mgxs_type)
mgxs.load_from_statepoint(statepoint)
mgxs.sparse = self.sparse
def get_mgxs(self, domain, mgxs_type):
"""Return the MGXS object for some domain and reaction rate type.

335
openmc/mgxs/mgxs.py
View file

@ -47,7 +47,7 @@ _DOMAINS = [openmc.Cell,
class MGXS(object):
"""An abstract multi-group cross section for some energy group structure
"""An abstract multi-group cross section for some energy group structure
within some spatial domain.
This class can be used for both OpenMC input generation and tally data
@ -89,13 +89,17 @@ class MGXS(object):
compute the cross section
tallies : OrderedDict
OpenMC tallies needed to compute the multi-group cross section
rxn_rate_tally : Tally
Derived tally for the reaction rate tally used in the numerator to
compute the multi-group cross section. This attribute is None
unless the multi-group cross section has been computed.
xs_tally : Tally
Derived tally for the multi-group cross section. This attribute
is None unless the multi-group cross section has been computed.
num_sumbdomains : Integral
num_subdomains : Integral
The number of subdomains is unity for 'material', 'cell' and 'universe'
domain types. When the This is equal to the number of cell instances
for 'distribcell' domain types (it is equal to unity prior to loading
for 'distribcell' domain types (it is equal to unity prior to loading
tally data from a statepoint file).
num_nuclides : Integral
The number of nuclides for which the multi-group cross section is
@ -103,6 +107,9 @@ class MGXS(object):
nuclides : list of str or 'sum'
A list of nuclide string names (e.g., 'U-238', 'O-16') when by_nuclide
is True and 'sum' when by_nuclide is False.
sparse : bool
Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format
for compressed data storage
"""
@ -120,7 +127,9 @@ class MGXS(object):
self._energy_groups = None
self._tally_trigger = None
self._tallies = None
self._rxn_rate_tally = None
self._xs_tally = None
self._sparse = False
self.name = name
self.by_nuclide = by_nuclide
@ -145,7 +154,9 @@ class MGXS(object):
clone._domain_type = self.domain_type
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
clone._xs_tally = copy.deepcopy(self.xs_tally, memo)
clone._rxn_rate_tally = copy.deepcopy(self._rxn_rate_tally, memo)
clone._xs_tally = copy.deepcopy(self._xs_tally, memo)
clone._sparse = self.sparse
clone._tallies = OrderedDict()
for tally_type, tally in self.tallies.items():
@ -195,10 +206,29 @@ class MGXS(object):
def tallies(self):
return self._tallies
@property
def rxn_rate_tally(self):
return self._rxn_rate_tally
@property
def xs_tally(self):
"""Computes multi-group cross section using OpenMC tally arithmetic."""
if self._xs_tally is None:
if self.tallies is None:
msg = 'Unable to get xs_tally since tallies have ' \
'not been loaded from a statepoint'
raise ValueError(msg)
self._xs_tally = self.rxn_rate_tally / self.tallies['flux']
self._compute_xs()
return self._xs_tally
@property
def sparse(self):
return self._sparse
@property
def num_subdomains(self):
tally = list(self.tallies.values())[0]
@ -249,6 +279,31 @@ class MGXS(object):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
self._tally_trigger = tally_trigger
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
# Sparsify or densify the derived MGXS tallies and the base tallies
if self._xs_tally:
self.xs_tally.sparse = sparse
if self._rxn_rate_tally:
self.rxn_rate_tally.sparse = sparse
for tally_name in self.tallies:
self.tallies[tally_name].sparse = sparse
self._sparse = sparse
@staticmethod
def get_mgxs(mgxs_type, domain=None, domain_type=None,
energy_groups=None, by_nuclide=False, name=''):
@ -466,8 +521,8 @@ class MGXS(object):
self.tallies[key].add_trigger(trigger_clone)
# Add all non-domain specific Filters (e.g., 'energy') to the Tally
for filter in filters:
self.tallies[key].add_filter(filter)
for add_filter in filters:
self.tallies[key].add_filter(add_filter)
# If this is a by-nuclide cross-section, add all nuclides to Tally
if self.by_nuclide and score != 'flux':
@ -477,8 +532,7 @@ class MGXS(object):
else:
self.tallies[key].add_nuclide('total')
@abc.abstractmethod
def compute_xs(self):
def _compute_xs(self):
"""Performs generic cleanup after a subclass' uses tally arithmetic to
compute a multi-group cross section as a derived tally.
@ -501,8 +555,9 @@ class MGXS(object):
self.xs_tally.add_nuclide(openmc.Nuclide(nuclide))
# Remove NaNs which may have resulted from divide-by-zero operations
self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean)
self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
self.xs_tally._mean = np.nan_to_num(self.xs_tally.mean)
self.xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
self.xs_tally.sparse = self.sparse
def load_from_statepoint(self, statepoint):
"""Extracts tallies in an OpenMC StatePoint with the data needed to
@ -565,11 +620,9 @@ class MGXS(object):
estimator=tally.estimator)
sp_tally = sp_tally.get_slice(tally.scores, filters,
filter_bins, tally.nuclides)
sp_tally.sparse = self.sparse
self.tallies[tally_type] = sp_tally
# Compute the cross section from the tallies
self.compute_xs()
def get_xs(self, groups='all', subdomains='all', nuclides='all',
xs_type='macro', order_groups='increasing', value='mean'):
"""Returns an array of multi-group cross sections.
@ -613,10 +666,6 @@ class MGXS(object):
"""
if self.xs_tally is None:
msg = 'Unable to get cross section since it has not been computed'
raise ValueError(msg)
cv.check_value('value', value, ['mean', 'std_dev', 'rel_err'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -701,11 +750,6 @@ class MGXS(object):
"""
if self.xs_tally is None:
msg = 'Unable to get a condensed coarse group cross section ' \
'since the fine group cross section has not been computed'
raise ValueError(msg)
cv.check_type('coarse_groups', coarse_groups, EnergyGroups)
cv.check_less_than('coarse groups', coarse_groups.num_groups,
self.num_groups, equality=True)
@ -716,7 +760,10 @@ class MGXS(object):
# Clone this MGXS to initialize the condensed version
condensed_xs = copy.deepcopy(self)
condensed_xs.energy_groups = coarse_groups
condensed_xs._rxn_rate_tally = None
condensed_xs._xs_tally = None
condensed_xs._sparse = False
condensed_xs._energy_groups = coarse_groups
# Build energy indices to sum across
energy_indices = []
@ -740,31 +787,29 @@ class MGXS(object):
std_dev = tally.get_reshaped_data(value='std_dev')
# Sum across all applicable fine energy group filters
for i, filter in enumerate(tally.filters):
if 'energy' not in filter.type:
for i, tally_filter in enumerate(tally.filters):
if 'energy' not in tally_filter.type:
continue
elif len(filter.bins) != len(fine_edges):
elif len(tally_filter.bins) != len(fine_edges):
continue
elif not np.allclose(filter.bins, fine_edges):
elif not np.allclose(tally_filter.bins, fine_edges):
continue
else:
filter.bins = coarse_groups.group_edges
tally_filter.bins = coarse_groups.group_edges
mean = np.add.reduceat(mean, energy_indices, axis=i)
std_dev = np.add.reduceat(std_dev**2, energy_indices, axis=i)
std_dev = np.sqrt(std_dev)
# Reshape condensed data arrays with one dimension for all filters
new_shape = \
(tally.num_filter_bins, tally.num_nuclides, tally.num_score_bins,)
mean = np.reshape(mean, new_shape)
std_dev = np.reshape(std_dev, new_shape)
mean = np.reshape(mean, tally.shape)
std_dev = np.reshape(std_dev, tally.shape)
# Override tally's data with the new condensed data
tally._mean = mean
tally._std_dev = std_dev
# Compute the energy condensed multi-group cross section
condensed_xs.compute_xs()
condensed_xs.sparse = self.sparse
return condensed_xs
def get_subdomain_avg_xs(self, subdomains='all'):
@ -792,11 +837,6 @@ class MGXS(object):
"""
if self.xs_tally is None:
msg = 'Unable to get subdomain-averaged cross section since the ' \
'subdomain-distributed cross section has not been computed'
raise ValueError(msg)
# Construct a collection of the subdomain filter bins to average across
if subdomains != 'all':
cv.check_iterable_type('subdomains', subdomains, Integral)
@ -807,8 +847,11 @@ class MGXS(object):
# Clone this MGXS to initialize the subdomain-averaged version
avg_xs = copy.deepcopy(self)
avg_xs._rxn_rate_tally = None
avg_xs._xs_tally = None
avg_xs._sparse = False
# If domain is distribcell, make the new domain 'cell'
# If domain is distribcell, make the new domain 'cell'
if self.domain_type == 'distribcell':
avg_xs.domain_type = 'cell'
@ -836,18 +879,15 @@ class MGXS(object):
domain_filter.num_bins = 1
# Reshape averaged data arrays with one dimension for all filters
new_shape = \
(tally.num_filter_bins, tally.num_nuclides, tally.num_score_bins,)
mean = np.reshape(mean, new_shape)
std_dev = np.reshape(std_dev, new_shape)
mean = np.reshape(mean, tally.shape)
std_dev = np.reshape(std_dev, tally.shape)
# Override tally's data with the new condensed data
tally._mean = mean
tally._std_dev = std_dev
# Compute the subdomain-averaged multi-group cross section
avg_xs.compute_xs()
avg_xs.sparse = self.sparse
return avg_xs
def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'):
@ -898,7 +938,7 @@ class MGXS(object):
string += '{0: <16}=\t{1}\n'.format('\tDomain ID', self.domain.id)
# If cross section data has not been computed, only print string header
if self.xs_tally is None:
if self.tallies is None:
print(string)
return
@ -985,11 +1025,6 @@ class MGXS(object):
"""
if self.xs_tally is None:
msg = 'Unable to get build HDF5 store since the ' \
'cross section has not been computed'
raise ValueError(msg)
import h5py
# Make directory if it does not exist
@ -1187,11 +1222,6 @@ class MGXS(object):
"""
if self.xs_tally is None:
msg = 'Unable to get Pandas DataFrame since the ' \
'cross section has not been computed'
raise ValueError(msg)
if groups != 'all':
cv.check_iterable_type('groups', groups, Integral)
if nuclides != 'all' and nuclides != 'sum':
@ -1312,13 +1342,12 @@ class TotalXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group total cross sections using OpenMC
tally arithmetic.
"""
self._xs_tally = self.tallies['total'] / self.tallies['flux']
super(TotalXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['total']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class TransportXS(MGXS):
@ -1359,18 +1388,21 @@ class TransportXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group transport cross sections using OpenMC
tally arithmetic."""
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
scatter_p1 = copy.deepcopy(self.tallies['scatter-P1'])
# Use tally slicing to remove scatter-P0 data from scatter-P1 tally
scatter_p1 = self.tallies['scatter-P1']
self.tallies['scatter-P1'] = scatter_p1.get_slice(scores=['scatter-P1'])
self.tallies['scatter-P1'].filters[-1].type = 'energy'
# Use tally slicing to remove scatter-P0 data from scatter-P1 tally
self.tallies['scatter-P1'] = \
scatter_p1.get_slice(scores=['scatter-P1'])
self._xs_tally = self.tallies['total'] - self.tallies['scatter-P1']
self._xs_tally /= self.tallies['flux']
super(TransportXS, self).compute_xs()
self.tallies['scatter-P1'].filters[-1].type = 'energy'
self._rxn_rate_tally = \
self.tallies['total'] - self.tallies['scatter-P1']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class AbsorptionXS(MGXS):
@ -1410,18 +1442,18 @@ class AbsorptionXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group absorption cross sections using OpenMC
tally arithmetic."""
self._xs_tally = self.tallies['absorption'] / self.tallies['flux']
super(AbsorptionXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['absorption']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class CaptureXS(MGXS):
"""A capture multi-group cross section.
The Neutron capture reaction rate is defined as the difference between
The neutron capture reaction rate is defined as the difference between
OpenMC's 'absorption' and 'fission' reaction rate score types. This includes
not only radiative capture, but all forms of neutron disappearance aside
from fission (e.g., MT > 100).
@ -1462,13 +1494,13 @@ class CaptureXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group capture cross sections using OpenMC
tally arithmetic."""
self._xs_tally = self.tallies['absorption'] - self.tallies['fission']
self._xs_tally /= self.tallies['flux']
super(CaptureXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = \
self.tallies['absorption'] - self.tallies['fission']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class FissionXS(MGXS):
@ -1508,12 +1540,12 @@ class FissionXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group fission cross sections using OpenMC
tally arithmetic."""
self._xs_tally = self.tallies['fission'] / self.tallies['flux']
super(FissionXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['fission']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class NuFissionXS(MGXS):
@ -1553,12 +1585,12 @@ class NuFissionXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the multi-group nu-fission cross sections using OpenMC
tally arithmetic."""
self._xs_tally = self.tallies['nu-fission'] / self.tallies['flux']
super(NuFissionXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['nu-fission']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class ScatterXS(MGXS):
@ -1598,12 +1630,12 @@ class ScatterXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the scattering multi-group cross sections using
OpenMC tally arithmetic."""
self._xs_tally = self.tallies['scatter'] / self.tallies['flux']
super(ScatterXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['scatter']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class NuScatterXS(MGXS):
@ -1643,12 +1675,12 @@ class NuScatterXS(MGXS):
return self._tallies
def compute_xs(self):
"""Computes the nu-scattering multi-group cross section using OpenMC
tally arithmetic."""
self._xs_tally = self.tallies['nu-scatter'] / self.tallies['flux']
super(NuScatterXS, self).compute_xs()
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['nu-scatter']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
class ScatterMatrixXS(MGXS):
@ -1710,29 +1742,30 @@ class ScatterMatrixXS(MGXS):
return self._tallies
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
# If using P0 correction subtract scatter-P1 from the diagonal
if self.correction == 'P0':
scatter_p1 = self.tallies['scatter-P1']
scatter_p1 = scatter_p1.get_slice(scores=['scatter-P1'])
energy_filter = self.tallies['scatter'].find_filter('energy')
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = self.tallies['scatter'] - scatter_p1
else:
self._rxn_rate_tally = self.tallies['scatter']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
@correction.setter
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))
self._correction = correction
def compute_xs(self):
"""Computes the multi-group scattering matrix using OpenMC
tally arithmetic."""
# If using P0 correction subtract scatter-P1 from the diagonal
if self.correction == 'P0':
scatter_p1 = self.tallies['scatter-P1']
scatter_p1 = scatter_p1.get_slice(scores=['scatter-P1'])
energy_filter = openmc.Filter(type='energy')
energy_filter.bins = self.energy_groups.group_edges
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
rxn_tally = self.tallies['scatter'] - scatter_p1
else:
rxn_tally = self.tallies['scatter']
self._xs_tally = rxn_tally / self.tallies['flux']
super(ScatterMatrixXS, self).compute_xs()
def get_xs(self, in_groups='all', out_groups='all',
subdomains='all', nuclides='all', xs_type='macro',
order_groups='increasing', value='mean'):
@ -1779,10 +1812,6 @@ class ScatterMatrixXS(MGXS):
"""
if self.xs_tally is None:
msg = 'Unable to get cross section since it has not been computed'
raise ValueError(msg)
cv.check_value('value', value, ['mean', 'std_dev', 'rel_err'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@ -1915,7 +1944,7 @@ class ScatterMatrixXS(MGXS):
string += '{0: <16}=\t{1}\n'.format('\tDomain ID', self.domain.id)
# If cross section data has not been computed, only print string header
if self.xs_tally is None:
if self.tallies is None:
print(string)
return
@ -2059,24 +2088,32 @@ class Chi(MGXS):
return self._tallies
def compute_xs(self):
@property
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
self._rxn_rate_tally = self.tallies['nu-fission-out']
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
@property
def xs_tally(self):
"""Computes chi fission spectrum using OpenMC tally arithmetic."""
# Retrieve the fission production tallies
nu_fission_in = self.tallies['nu-fission-in']
nu_fission_out = self.tallies['nu-fission-out']
if self._xs_tally is None:
nu_fission_in = self.tallies['nu-fission-in']
# Remove coarse energy filter to keep it out of tally arithmetic
energy_filter = nu_fission_in.find_filter('energy')
nu_fission_in.remove_filter(energy_filter)
# Remove coarse energy filter to keep it out of tally arithmetic
energy_filter = nu_fission_in.find_filter('energy')
nu_fission_in.remove_filter(energy_filter)
# Compute chi
self._xs_tally = nu_fission_out / nu_fission_in
# Compute chi
self._xs_tally = self.rxn_rate_tally / nu_fission_in
super(Chi, self)._compute_xs()
# Add the coarse energy filter back to the nu-fission tally
nu_fission_in.add_filter(energy_filter)
# Add the coarse energy filter back to the nu-fission tally
nu_fission_in.add_filter(energy_filter)
super(Chi, self).compute_xs()
return self._xs_tally
def get_xs(self, groups='all', subdomains='all', nuclides='all',
xs_type='macro', order_groups='increasing', value='mean'):
@ -2121,10 +2158,6 @@ class Chi(MGXS):
"""
if self.xs_tally is None:
msg = 'Unable to get cross section since it has not been computed'
raise ValueError(msg)
cv.check_value('value', value, ['mean', 'std_dev', 'rel_err'])
cv.check_value('xs_type', xs_type, ['macro', 'micro'])

28
openmc/opencg_compatible.py
View file

@ -725,11 +725,11 @@ def get_openmc_cell(opencg_cell):
else:
openmc_cell.fill = get_openmc_material(fill)
if opencg_cell.rotation:
if opencg_cell.rotation is not None:
rotation = np.asarray(opencg_cell.rotation, dtype=np.float64)
openmc_cell.rotation = rotation
if opencg_cell.translation:
if opencg_cell.translation is not None:
translation = np.asarray(opencg_cell.translation, dtype=np.float64)
openmc_cell.translation = translation
@ -881,16 +881,30 @@ def get_opencg_lattice(openmc_lattice):
universes = openmc_lattice.universes
outer = openmc_lattice.outer
# Convert 2D dimension to 3D for OpenCG
if len(dimension) == 2:
new_dimension = np.ones(3, dtype=np.int)
new_dimension[:2] = dimension
dimension = new_dimension
# Convert 2D pitch to 3D for OpenCG
if len(pitch) == 2:
new_pitch = np.ones(3, dtype=np.float64) * np.inf
new_pitch = np.ones(3, dtype=np.float64) * np.finfo(np.float64).max
new_pitch[:2] = pitch
pitch = new_pitch
# Convert 2D lower left to 3D for OpenCG
if len(lower_left) == 2:
new_lower_left = np.ones(3, dtype=np.float64)
new_lower_left = np.ones(3, dtype=np.float64) * np.finfo(np.float64).min
new_lower_left[:2] = lower_left
lower_left = new_lower_left
# Convert 2D universes array to 3D for OpenCG
if len(universes.shape) == 2:
new_universes = universes.copy()
new_universes.shape = (1,) + universes.shape
universes = new_universes
# Initialize an empty array for the OpenCG nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
dtype=opencg.Universe)
@ -905,7 +919,7 @@ def get_opencg_lattice(openmc_lattice):
for z in range(dimension[2]):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[x][dimension[1]-y-1][z].id
universe_id = universes[z][y][x].id
universe_array[z][y][x] = unique_universes[universe_id]
opencg_lattice = opencg.Lattice(lattice_id, name)
@ -963,7 +977,7 @@ def get_openmc_lattice(opencg_lattice):
outer = opencg_lattice.outside
# Initialize an empty array for the OpenMC nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)),
universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
dtype=openmc.Universe)
# Create OpenMC Universes for each unique nested Universe in this Lattice
@ -977,7 +991,7 @@ def get_openmc_lattice(opencg_lattice):
for y in range(dimension[1]):
for x in range(dimension[0]):
universe_id = universes[z][y][x].id
universe_array[x][y][z] = unique_universes[universe_id]
universe_array[z][y][x] = unique_universes[universe_id]
# Reverse y-dimension in array to match ordering in OpenCG
universe_array = universe_array[:, ::-1, :]

212
openmc/settings.py
View file

@ -9,6 +9,7 @@ import numpy as np
from openmc.clean_xml import *
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
from openmc.source import Source
if sys.version_info[0] >= 3:
basestring = str
@ -36,8 +37,8 @@ class SettingsFile(object):
type are 'variance', 'std_dev', and 'rel_err'. The threshold value
should be a float indicating the variance, standard deviation, or
relative error used.
source_file : str
Path to a source file
source : Iterable of openmc.source.Source
Distribution of source sites in space, angle, and energy
output : dict
Dictionary indicating what files to output. Valid keys are 'summary',
'cross_sections', 'tallies', and 'distribmats'. Values corresponding to
@ -134,14 +135,7 @@ class SettingsFile(object):
self._keff_trigger = None
# Source subelement
self._source_subelement = None
self._source_file = None
self._source_space_type = None
self._source_space_params = None
self._source_angle_type = None
self._source_angle_params = None
self._source_energy_type = None
self._source_energy_params = None
self._source = None
self._confidence_intervals = None
self._cross_sections = None
@ -228,32 +222,8 @@ class SettingsFile(object):
return self._keff_trigger
@property
def source_file(self):
return self._source_file
@property
def source_space_type(self):
return self._source_space_type
@property
def source_space_params(self):
return self._source_space_params
@property
def source_angle_type(self):
return self._source_angle_type
@property
def source_angle_params(self):
return self._source_angle_params
@property
def source_energy_type(self):
return self._source_energy_type
@property
def source_energy_params(self):
return self._source_energy_params
def source(self):
return self._source
@property
def confidence_intervals(self):
@ -468,124 +438,13 @@ class SettingsFile(object):
self._keff_trigger = keff_trigger
@source_file.setter
def source_file(self, source_file):
check_type('source file', source_file, basestring)
self._source_file = source_file
def set_source_space(self, stype, params):
"""Defined the spatial bounds of the external/starting source.
Parameters
----------
stype : str
The type of spatial distribution. Valid options are "box",
"fission", and "point". A "box" spatial distribution has coordinates
sampled uniformly in a parallelepiped. A "fission" spatial
distribution samples locations from a "box" distribution but only
locations in fissionable materials are accepted. A "point" spatial
distribution has coordinates specified by a triplet.
params : Iterable of float
For a "box" or "fission" spatial distribution, ``params`` should be
given as six real numbers, the first three of which specify the
lower-left corner of a parallelepiped and the last three of which
specify the upper-right corner. Source sites are sampled uniformly
through that parallelepiped.
For a "point" spatial distribution, ``params`` should be given as
three real numbers which specify the (x,y,z) location of an
isotropic point source
"""
check_type('source space type', stype, basestring)
check_value('source space type', stype, ['box', 'fission', 'point'])
check_type('source space parameters', params, Iterable, Real)
if stype in ['box', 'fission']:
check_length('source space parameters for a '
'box/fission distribution', params, 6)
elif stype == 'point':
check_length('source space parameters for a point source',
params, 3)
self._source_space_type = stype
self._source_space_params = params
def set_source_angle(self, stype, params=[]):
"""Defined the angular distribution of the external/starting source.
Parameters
----------
stype : str
The type of angular distribution. Valid options are "isotropic" and
"monodirectional". The angle of the particle emitted from a source
site is isotropic if the "isotropic" option is given. The angle of
the particle emitted from a source site is the direction specified
in ``params`` if the "monodirectional" option is given.
params : Iterable of float
For an "isotropic" angular distribution, ``params`` should not
be specified.
For a "monodirectional" angular distribution, ``params`` should
be given as three floats which specify the angular cosines
with respect to each axis.
"""
check_type('source angle type', stype, basestring)
check_value('source angle type', stype,
['isotropic', 'monodirectional'])
check_type('source angle parameters', params, Iterable, Real)
if stype == 'isotropic' and params is not None:
msg = 'Unable to set source angle parameters since they are not ' \
'it is not supported for isotropic type sources'
raise ValueError(msg)
elif stype == 'monodirectional':
check_length('source angle parameters for a monodirectional '
'source', params, 3)
self._source_angle_type = stype
self._source_angle_params = params
def set_source_energy(self, stype, params=[]):
"""Defined the energy distribution of the external/starting source.
Parameters
----------
stype : str
The type of energy distribution. Valid options are "monoenergetic",
"watt", and "maxwell". The "monoenergetic" option produces source
sites at a single energy. The "watt" option produces source sites
whose energy is sampled from a Watt fission spectrum. The "maxwell"
option produce source sites whose energy is sampled from a Maxwell
fission spectrum.
params : Iterable of float
For a "monoenergetic" energy distribution, ``params`` should be
given as the energy in MeV of the source sites.
For a "watt" energy distribution, ``params`` should be given as two
real numbers :math:`a` and :math:`b` that parameterize the
distribution :math:`p(E) dE = c e^{-E/a} \sinh \sqrt{b \, E} dE`.
For a "maxwell" energy distribution, ``params`` should be given as
one real number :math:`a` that parameterizes the distribution
:math:`p(E) dE = c E e^{-E/a} dE`.
"""
check_type('source energy type', stype, basestring)
check_value('source energy type', stype,
['monoenergetic', 'watt', 'maxwell'])
check_type('source energy parameters', params, Iterable, Real)
if stype in ['monoenergetic', 'maxwell']:
check_length('source energy parameters for a monoenergetic '
'or Maxwell source', params, 1)
elif stype == 'watt':
check_length('source energy parameters for a Watt source',
params, 2)
self._source_energy_type = stype
self._source_energy_params = params
@source.setter
def source(self, source):
if isinstance(source, Source):
self._source = [source,]
else:
check_type('source distribution', source, Iterable, Source)
self._source = source
@output.setter
def output(self, output):
@ -924,45 +783,9 @@ class SettingsFile(object):
subelement.text = str(self._keff_trigger[key]).lower()
def _create_source_subelement(self):
self._create_source_space_subelement()
self._create_source_energy_subelement()
self._create_source_angle_subelement()
def _create_source_space_subelement(self):
if self._source_space_params is not None:
if self._source_subelement is None:
self._source_subelement = ET.SubElement(self._settings_file,
"source")
element = ET.SubElement(self._source_subelement, "space")
element.set("type", self._source_space_type)
subelement = ET.SubElement(element, "parameters")
subelement.text = ' '.join(map(str, self._source_space_params))
def _create_source_angle_subelement(self):
if self._source_angle_params is not None:
if self._source_subelement is None:
self._source_subelement = ET.SubElement(self._settings_file,
"source")
element = ET.SubElement(self._source_subelement, "angle")
element.set("type", self._source_angle_type)
subelement = ET.SubElement(element, "parameters")
subelement.text = ' '.join(map(str, self._source_angle_params))
def _create_source_energy_subelement(self):
if self._source_energy_params is not None:
if self._source_subelement is None:
self._source_subelement = ET.SubElement(self._settings_file,
"source")
element = ET.SubElement(self._source_subelement, "energy")
element.set("type", self._source_energy_type)
subelement = ET.SubElement(element, "parameters")
subelement.text = ' '.join(map(str, self._source_energy_params))
if self.source is not None:
for source in self.source:
self._settings_file.append(source.to_xml())
def _create_output_subelement(self):
if self._output is not None:
@ -1012,16 +835,19 @@ class SettingsFile(object):
# Separate subelement
if self._sourcepoint_separate is not None:
element = ET.SubElement(self._settings_file, "source_point")
subelement = ET.SubElement(element, "separate")
subelement.text = str(self._sourcepoint_separate).lower()
# Write subelement
if self._sourcepoint_write is not None:
element = ET.SubElement(self._settings_file, "source_point")
subelement = ET.SubElement(element, "write")
subelement.text = str(self._sourcepoint_write).lower()
# Overwrite latest subelement
if self._sourcepoint_overwrite is not None:
element = ET.SubElement(self._settings_file, "source_point")
subelement = ET.SubElement(element, "overwrite_latest")
subelement.text = str(self._sourcepoint_overwrite).lower()

117
openmc/source.py Normal file
View file

@ -0,0 +1,117 @@
from numbers import Real
import sys
from xml.etree import ElementTree as ET
from openmc.stats.univariate import Univariate
from openmc.stats.multivariate import UnitSphere, Spatial
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
class Source(object):
"""Distribution of phase space coordinates for source sites.
Parameters
----------
space : openmc.stats.Spatial, optional
Spatial distribution of source sites
angle : openmc.stats.UnitSphere, optional
Angular distribution of source sites
energy : openmc.stats.Univariate, optional
Energy distribution of source sites
filename : str, optional
Source file from which sites should be sampled
strength : Real
Strength of the source
Attributes
----------
space : openmc.stats.Spatial or None
Spatial distribution of source sites
angle : openmc.stats.UnitSphere or None
Angular distribution of source sites
energy : openmc.stats.Univariate or None
Energy distribution of source sites
file : str or None
Source file from which sites should be sampled
strength : Real
Strength of the source
"""
def __init__(self, space=None, angle=None, energy=None, filename=None, strength=1.0):
self._space = None
self._angle = None
self._energy = None
self._file = None
if space is not None:
self.space = space
if angle is not None:
self.angle = angle
if energy is not None:
self.energy = energy
if filename is not None:
self.file = filename
self.strength = strength
@property
def file(self):
return self._file
@property
def space(self):
return self._space
@property
def angle(self):
return self._angle
@property
def energy(self):
return self._energy
@property
def strength(self):
return self._strength
@file.setter
def file(self, filename):
cv.check_type('source file', filename, basestring)
self._file = filename
@space.setter
def space(self, space):
cv.check_type('spatial distribution', space, Spatial)
self._space = space
@angle.setter
def angle(self, angle):
cv.check_type('angular distribution', angle, UnitSphere)
self._angle = angle
@energy.setter
def energy(self, energy):
cv.check_type('energy distribution', energy, Univariate)
self._energy = energy
@strength.setter
def strength(self, strength):
cv.check_type('source strength', strength, Real)
cv.check_greater_than('source strength', strength, 0.0, True)
self._strength = strength
def to_xml(self):
element = ET.Element("source")
element.set("strength", str(self.strength))
if self.file is not None:
element.set("file", self.file)
if self.space is not None:
element.append(self.space.to_xml())
if self.angle is not None:
element.append(self.angle.to_xml())
if self.energy is not None:
element.append(self.energy.to_xml('energy'))
return element

109
openmc/statepoint.py
View file

@ -3,15 +3,16 @@ import re
import numpy as np
import openmc
import openmc.checkvalue as cv
if sys.version > '3':
long = int
class StatePoint(object):
"""State information on a simulation at a certain point in time (at the end of a
given batch). Statepoints can be used to analyze tally results as well as
restart a simulation.
"""State information on a simulation at a certain point in time (at the end
of a given batch). Statepoints can be used to analyze tally results as well
as restart a simulation.
Attributes
----------
@ -75,6 +76,9 @@ class StatePoint(object):
energy of the source site.
source_present : bool
Indicate whether source sites are present
sparse : bool
Whether or not the tallies uses SciPy's LIL sparse matrix format for
compressed data storage
tallies : dict
Dictionary whose keys are tally IDs and whose values are Tally objects
tallies_present : bool
@ -110,6 +114,7 @@ class StatePoint(object):
self._tallies_read = False
self._summary = False
self._global_tallies = None
self._sparse = False
def close(self):
self._f.close()
@ -318,6 +323,10 @@ class StatePoint(object):
def source_present(self):
return self._f['source_present'].value > 0
@property
def sparse(self):
return self._sparse
@property
def tallies(self):
if not self._tallies_read:
@ -340,7 +349,8 @@ class StatePoint(object):
for tally_key in tally_keys:
# Read the Tally size specifications
n_realizations = self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
n_realizations = \
self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id=tally_key)
@ -373,7 +383,8 @@ class StatePoint(object):
'variable')
# Read the number of Filters
n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value
n_filters = \
self._f['{0}{1}/n_filters'.format(base, tally_key)].value
subbase = '{0}{1}/filter '.format(base, tally_key)
@ -381,10 +392,8 @@ class StatePoint(object):
for j in range(1, n_filters+1):
# Read the Filter type
filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode()
# Read the Filter offset
offset = self._f['{0}{1}/offset'.format(subbase, j)].value
filter_type = \
self._f['{0}{1}/type'.format(subbase, j)].value.decode()
n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value
@ -392,47 +401,42 @@ class StatePoint(object):
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
# Create Filter object
filter = openmc.Filter(filter_type, bins)
filter.offset = offset
filter.num_bins = n_bins
new_filter = openmc.Filter(filter_type, bins)
new_filter.num_bins = n_bins
if filter_type == 'mesh':
mesh_ids = self._f['tallies/meshes/ids'].value
mesh_keys = self._f['tallies/meshes/keys'].value
key = mesh_keys[mesh_ids == bins][0]
filter.mesh = self.meshes[key]
new_filter.mesh = self.meshes[key]
# Add Filter to the Tally
tally.add_filter(filter)
tally.add_filter(new_filter)
# Read Nuclide bins
nuclide_names = self._f['{0}{1}/nuclides'.format(base, tally_key)].value
nuclide_names = \
self._f['{0}{1}/nuclides'.format(base, tally_key)].value
# Add all Nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.add_nuclide(nuclide)
# Read score bins
n_score_bins = self._f['{0}{1}/n_score_bins'.format(base, tally_key)].value
tally.num_score_bins = n_score_bins
scores = self._f['{0}{1}/score_bins'.format(
base, tally_key)].value
n_user_scores = self._f['{0}{1}/n_user_score_bins'
.format(base, tally_key)].value
n_score_bins = self._f['{0}{1}/n_score_bins'
.format(base, tally_key)].value
# Compute and set the filter strides
for i in range(n_filters):
filter = tally.filters[i]
filter.stride = n_score_bins * len(nuclide_names)
tally_filter = tally.filters[i]
tally_filter.stride = n_score_bins * len(nuclide_names)
for j in range(i+1, n_filters):
filter.stride *= tally.filters[j].num_bins
tally_filter.stride *= tally.filters[j].num_bins
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = self._f['{0}{1}/moment_orders'.format(
base, tally_key)].value
@ -447,6 +451,7 @@ class StatePoint(object):
tally.add_score(score)
# Add Tally to the global dictionary of all Tallies
tally.sparse = self.sparse
self._tallies[tally_key] = tally
self._tallies_read = True
@ -471,6 +476,26 @@ class StatePoint(object):
def with_summary(self):
return False if self.summary is None else True
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within each Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
self._sparse = sparse
# Update tally sparsities
if self._tallies_read:
for tally_id in self.tallies:
self.tallies[tally_id].sparse = self.sparse
def get_tally(self, scores=[], filters=[], nuclides=[],
name=None, id=None, estimator=None):
"""Finds and returns a Tally object with certain properties.
@ -543,13 +568,13 @@ class StatePoint(object):
contains_filters = True
# Iterate over the Filters requested by the user
for filter in filters:
for outer_filter in filters:
contains_filters = False
# Test if requested filter is a subset of any of the test
# tally's filters and if so continue to next filter
for test_filter in test_tally.filters:
if test_filter.is_subset(filter):
for inner_filter in test_tally.filters:
if inner_filter.is_subset(outer_filter):
contains_filters = True
break
@ -615,29 +640,29 @@ class StatePoint(object):
tally.name = summary.tallies[tally_id].name
tally.with_summary = True
for filter in tally.filters:
if filter.type == 'surface':
for tally_filter in tally.filters:
if tally_filter.type == 'surface':
surface_ids = []
for bin in filter.bins:
for bin in tally_filter.bins:
surface_ids.append(summary.surfaces[bin].id)
filter.bins = surface_ids
tally_filter.bins = surface_ids
if filter.type in ['cell', 'distribcell']:
if tally_filter.type in ['cell', 'distribcell']:
distribcell_ids = []
for bin in filter.bins:
for bin in tally_filter.bins:
distribcell_ids.append(summary.cells[bin].id)
filter.bins = distribcell_ids
tally_filter.bins = distribcell_ids
if filter.type == 'universe':
if tally_filter.type == 'universe':
universe_ids = []
for bin in filter.bins:
for bin in tally_filter.bins:
universe_ids.append(summary.universes[bin].id)
filter.bins = universe_ids
tally_filter.bins = universe_ids
if filter.type == 'material':
if tally_filter.type == 'material':
material_ids = []
for bin in filter.bins:
for bin in tally_filter.bins:
material_ids.append(summary.materials[bin].id)
filter.bins = material_ids
tally_filter.bins = material_ids
self._summary = summary

2
openmc/stats/__init__.py Normal file
View file

@ -0,0 +1,2 @@
from openmc.stats.univariate import *
from openmc.stats.multivariate import *

360
openmc/stats/multivariate.py Normal file
View file

@ -0,0 +1,360 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from math import pi
from numbers import Real
import sys
from xml.etree import ElementTree as ET
import numpy as np
import openmc.checkvalue as cv
from openmc.stats.univariate import Univariate, Uniform
if sys.version_info[0] >= 3:
basestring = str
class UnitSphere(object):
"""Distribution of points on the unit sphere.
This abstract class is used for angular distributions, since a direction is
represented as a unit vector (i.e., vector on the unit sphere).
Parameters
----------
reference_uvw : Iterable of Real
Direction from which polar angle is measured
Attributes
----------
reference_uvw : Iterable of Real
Direction from which polar angle is measured
"""
__metaclass__ = ABCMeta
def __init__(self, reference_uvw=None):
self._reference_uvw = None
if reference_uvw is not None:
self.reference_uvw = reference_uvw
@property
def reference_uvw(self):
return self._reference_uvw
@reference_uvw.setter
def reference_uvw(self, uvw):
cv.check_type('reference direction', uvw, Iterable, Real)
uvw = np.asarray(uvw)
self._reference_uvw = uvw/np.linalg.norm(uvw)
@abstractmethod
def to_xml(self):
return ''
class PolarAzimuthal(UnitSphere):
"""Angular distribution represented by polar and azimuthal angles
This distribution allows one to specify the distribution of the cosine of
the polar angle and the azimuthal angle independently of once another.
Parameters
----------
mu : Univariate
Distribution of the cosine of the polar angle
phi : Univariate
Distribution of the azimuthal angle in radians
reference_uvw : Iterable of Real
Direction from which polar angle is measured. Defaults to the positive
z-direction.
Attributes
----------
mu : Univariate
Distribution of the cosine of the polar angle
phi : Univariate
Distribution of the azimuthal angle in radians
"""
def __init__(self, mu=None, phi=None, reference_uvw=[0., 0., 1.]):
super(PolarAzimuthal, self).__init__(reference_uvw)
if mu is not None:
self.mu = mu
else:
self.mu = Uniform(-1., 1.)
if phi is not None:
self.phi = phi
else:
self.phi = Uniform(0., 2*pi)
@property
def mu(self):
return self._mu
@property
def phi(self):
return self._phi
@mu.setter
def mu(self, mu):
cv.check_type('cosine of polar angle', mu, Univariate)
self._mu = mu
@phi.setter
def phi(self, phi):
cv.check_type('azimuthal angle', phi, Univariate)
self._phi = phi
def to_xml(self):
element = ET.Element('angle')
element.set("type", "mu-phi")
if self.reference_uvw is not None:
element.set("reference_uvw", ' '.join(map(str, self.reference_uvw)))
element.append(self.mu.to_xml('mu'))
element.append(self.phi.to_xml('phi'))
return element
class Isotropic(UnitSphere):
"""Isotropic angular distribution.
"""
def __init__(self):
super(Isotropic, self).__init__()
def to_xml(self):
element = ET.Element('angle')
element.set("type", "isotropic")
return element
class Monodirectional(UnitSphere):
"""Monodirectional angular distribution.
A monodirectional angular distribution is one for which the polar and
azimuthal angles are always the same. It is completely specified by the
reference direction vector.
Parameters
----------
reference_uvw : Iterable of Real
Direction from which polar angle is measured. Defaults to the positive
x-direction.
"""
def __init__(self, reference_uvw=[1., 0., 0.]):
super(Monodirectional, self).__init__(reference_uvw)
def to_xml(self):
element = ET.Element('angle')
element.set("type", "monodirectional")
if self.reference_uvw is not None:
element.set("reference_uvw", ' '.join(map(str, self.reference_uvw)))
return element
class Spatial(object):
"""Distribution of locations in three-dimensional Euclidean space.
Classes derived from this abstract class can be used for spatial
distributions of source sites.
"""
__metaclass__ = ABCMeta
def __init__(self):
pass
@abstractmethod
def to_xml(self):
return ''
class CartesianIndependent(Spatial):
"""Spatial distribution with independent x, y, and z distributions.
This distribution allows one to specify a coordinates whose x-, y-, and z-
components are sampled independently from one another.
Parameters
----------
x : Univariate
Distribution of x-coordinates
y : Univariate
Distribution of y-coordinates
z : Univariate
Distribution of z-coordinates
Attributes
----------
x : Univariate
Distribution of x-coordinates
y : Univariate
Distribution of y-coordinates
z : Univariate
Distribution of z-coordinates
"""
def __init__(self, x, y, z):
super(CartesianIndependent, self).__init__()
self.x = x
self.y = y
self.z = z
@property
def x(self):
return self._x
@property
def y(self):
return self._y
@property
def z(self):
return self._z
@x.setter
def x(self, x):
cv.check_type('x coordinate', x, Univariate)
self._x = x
@y.setter
def y(self, y):
cv.check_type('y coordinate', y, Univariate)
self._y = y
@z.setter
def z(self, z):
cv.check_type('z coordinate', z, Univariate)
self._z = z
def to_xml(self):
element = ET.Element('space')
element.set('type', 'cartesian')
element.append(self.x.to_xml('x'))
element.append(self.y.to_xml('y'))
element.append(self.z.to_xml('z'))
return element
class Box(Spatial):
"""Uniform distribution of coordinates in a rectangular cuboid.
Parameters
----------
lower_left : Iterable of Real
Lower-left coordinates of cuboid
upper_right : Iterable of Real
Upper-right coordinates of cuboid
only_fissionable : bool, optional
Whether spatial sites should only be accepted if they occur in
fissionable materials
Attributes
----------
lower_left : Iterable of Real
Lower-left coordinates of cuboid
upper_right : Iterable of Real
Upper-right coordinates of cuboid
only_fissionable : bool, optional
Whether spatial sites should only be accepted if they occur in
fissionable materials
"""
def __init__(self, lower_left, upper_right, only_fissionable=False):
super(Box, self).__init__()
self.lower_left = lower_left
self.upper_right = upper_right
self.only_fissionable = only_fissionable
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
return self._upper_right
@property
def only_fissionable(self):
return self._only_fissionable
@lower_left.setter
def lower_left(self, lower_left):
cv.check_type('lower left coordinate', lower_left, Iterable, Real)
cv.check_length('lower left coordinate', lower_left, 3)
self._lower_left = lower_left
@upper_right.setter
def upper_right(self, upper_right):
cv.check_type('upper right coordinate', upper_right, Iterable, Real)
cv.check_length('upper right coordinate', upper_right, 3)
self._upper_right = upper_right
@only_fissionable.setter
def only_fissionable(self, only_fissionable):
cv.check_type('only fissionable', only_fissionable, bool)
self._only_fissionable = only_fissionable
def to_xml(self):
element = ET.Element('space')
if self.only_fissionable:
element.set("type", "fission")
else:
element.set("type", "box")
params = ET.SubElement(element, "parameters")
params.text = ' '.join(map(str, self.lower_left)) + ' ' + \
' '.join(map(str, self.upper_right))
return element
class Point(Spatial):
"""Delta function in three dimensions.
This spatial distribution can be used for a point source where sites are
emitted at a specific location given by its Cartesian coordinates.
Parameters
----------
xyz : Iterable of Real
Cartesian coordinates of location
Attributes
----------
xyz : Iterable of Real
Cartesian coordinates of location
"""
def __init__(self, xyz):
super(Point, self).__init__()
self.xyz = xyz
@property
def xyz(self):
return self._xyz
@xyz.setter
def xyz(self, xyz):
cv.check_type('coordinate', xyz, Iterable, Real)
cv.check_length('coordinate', xyz, 3)
self._xyz = xyz
def to_xml(self):
element = ET.Element('space')
element.set("type", "point")
params = ET.SubElement(element, "parameters")
params.text = ' '.join(map(str, self.xyz))
return element

310
openmc/stats/univariate.py Normal file
View file

@ -0,0 +1,310 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
from numbers import Real
import sys
from xml.etree import ElementTree as ET
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
class Univariate(object):
"""Probability distribution of a single random variable.
The Univariate class is an abstract class that can be derived to implement a
specific probability distribution.
"""
__metaclass__ = ABCMeta
def __init__(self):
pass
@abstractmethod
def to_xml(self):
return ''
class Discrete(Univariate):
"""Distribution characterized by a probability mass function.
The Discrete distribution assigns probability values to discrete values of a
random variable, rather than expressing the distribution as a continuous
random variable.
Parameters
----------
x : Iterable of Real
Values of the random variable
p : Iterable of Real
Discrete probability for each value
Attributes
----------
x : Iterable of Real
Values of the random variable
p : Iterable of Real
Discrete probability for each value
"""
def __init__(self, x, p):
super(Discrete, self).__init__()
self.x = x
self.p = p
@property
def x(self):
return self._x
@property
def p(self):
return self._p
@x.setter
def x(self, x):
if cv._isinstance(x, Real):
x = [x]
cv.check_type('discrete values', x, Iterable, Real)
self._x = x
@p.setter
def p(self, p):
if cv._isinstance(p, Real):
p = [p]
cv.check_type('discrete probabilities', p, Iterable, Real)
for pk in p:
cv.check_greater_than('discrete probability', pk, 0.0, True)
self._p = p
def to_xml(self, element_name):
element = ET.Element(element_name)
element.set("type", "discrete")
params = ET.SubElement(element, "parameters")
params.text = ' '.join(map(str, self.x)) + ' ' + ' '.join(map(str, self.p))
return element
class Uniform(Univariate):
"""Distribution with constant probability over a finite interval [a,b]
Parameters
----------
a : float, optional
Lower bound of the sampling interval. Defaults to zero.
b : float, optional
Upper bound of the sampling interval. Defaults to unity.
Attributes
----------
a : float
Lower bound of the sampling interval
b : float
Upper bound of the sampling interval
"""
def __init__(self, a=0.0, b=1.0):
super(Uniform, self).__init__()
self.a = a
self.b = b
@property
def a(self):
return self._a
@property
def b(self):
return self._b
@a.setter
def a(self, a):
cv.check_type('Uniform a', a, Real)
self._a = a
@b.setter
def b(self, b):
cv.check_type('Uniform b', b, Real)
self._b = b
def to_xml(self, element_name):
element = ET.Element(element_name)
element.set("type", "uniform")
element.set("parameters", '{} {}'.format(self.a, self.b))
return element
class Maxwell(Univariate):
"""Maxwellian distribution in energy.
The Maxwellian distribution in energy is characterized by a single parameter
:math:`\theta` and has a density function :math:`p(E) dE = c E e^{-E/\theta}
dE`.
Parameters
----------
theta : float
Effective temperature for distribution
Attributes
----------
theta : float
Effective temperature for distribution
"""
def __init__(self, theta):
super(Maxwell, self).__init__()
self.theta = theta
@property
def theta(self):
return self._theta
@theta.setter
def theta(self, theta):
cv.check_type('Maxwell temperature', theta, Real)
cv.check_greater_than('Maxwell temperature', theta, 0.0)
self._theta = theta
def to_xml(self, element_name):
element = ET.Element(element_name)
element.set("type", "maxwell")
element.set("parameters", str(self.theta))
return element
class Watt(Univariate):
"""Watt fission energy spectrum.
The Watt fission energy spectrum is characterized by two parameters
:math:`a` and :math:`b` and has density function :math:`p(E) dE = c e^{-E/a}
\sinh \sqrt{b \, E} dE`.
Parameters
----------
a : float
First parameter of distribution
b : float
Second parameter of distribution
Attributes
----------
a : float
First parameter of distribution
b : float
Second parameter of distribution
"""
def __init__(self, a=0.988, b=2.249):
super(Watt, self).__init__()
self.a = a
self.b = b
@property
def a(self):
return self._a
@property
def b(self):
return self._b
@a.setter
def a(self, a):
cv.check_type('Watt a', a, Real)
cv.check_greater_than('Watt a', a, 0.0)
self._a = a
@b.setter
def b(self, b):
cv.check_type('Watt b', b, Real)
cv.check_greater_than('Watt b', b, 0.0)
self._b = b
def to_xml(self, element_name):
element = ET.Element(element_name)
element.set("type", "watt")
element.set("parameters", '{} {}'.format(self.a, self.b))
return element
class Tabular(Univariate):
"""Piecewise continuous probability distribution.
This class is used to represent a probability distribution whose density
function is tabulated at specific values and is either histogram or linearly
interpolated between points.
Parameters
----------
x : Iterable of Real
Tabulated values of the random variable
p : Iterable of Real
Tabulated probabilities
interpolation : {'histogram', 'linear-linear'}, optional
Indicate whether the density function is constant between tabulated
points or linearly-interpolated.
Attributes
----------
x : Iterable of Real
Tabulated values of the random variable
p : Iterable of Real
Tabulated probabilities
interpolation : {'histogram', 'linear-linear'}, optional
Indicate whether the density function is constant between tabulated
points or linearly-interpolated.
"""
def __init__(self, x, p, interpolation='linear-linear'):
super(Tabular, self).__init__()
self.x = x
self.p = p
self.interpolation = interpolation
@property
def x(self):
return self._x
@property
def p(self):
return self._p
@property
def interpolation(self):
return self._interpolation
@x.setter
def x(self, x):
cv.check_type('tabulated values', x, Iterable, Real)
self._x = x
@p.setter
def p(self, p):
cv.check_type('tabulated probabilities', p, Iterable, Real)
for pk in p:
cv.check_greater_than('tabulated probability', pk, 0.0, True)
self._p = p
@interpolation.setter
def interpolation(self, interpolation):
cv.check_value('interpolation', interpolation,
['linear-linear', 'histogram'])
self._interpolation = interpolation
def to_xml(self, element_name):
element = ET.Element(element_name)
element.set("type", "tabular")
element.set("interpolation", self.interpolation)
params = ET.SubElement(element, "parameters")
params.text = ' '.join(map(str, self.x)) + ' ' + ' '.join(map(str, self.p))
return element

64
openmc/summary.py
View file

@ -1,4 +1,6 @@
from collections import Iterable
import numpy as np
import re
import openmc
from openmc.region import Region
@ -271,6 +273,11 @@ class Summary(object):
cell.region = Region.from_expression(
region, {s.id: s for s in self.surfaces.values()})
# Get the distribcell index
ind = self._f['geometry/cells'][key]['distribcell_index'].value
if ind != 0:
cell.distribcell_index = ind
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell
@ -329,11 +336,8 @@ class Summary(object):
self._f['geometry/lattices'][key]['lower_left'][...]
pitch = self._f['geometry/lattices'][key]['pitch'][...]
outer = self._f['geometry/lattices'][key]['outer'].value
universe_ids = \
self._f['geometry/lattices'][key]['universes'][...]
universe_ids = np.swapaxes(universe_ids, 0, 1)
universe_ids = np.swapaxes(universe_ids, 1, 2)
self._f['geometry/lattices'][key]['universes'][...]
# Create the Lattice
lattice = openmc.RectLattice(lattice_id=lattice_id, name=name)
@ -349,22 +353,22 @@ class Summary(object):
universes = \
np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
for x in range(universe_ids.shape[0]):
for z in range(universe_ids.shape[0]):
for y in range(universe_ids.shape[1]):
for z in range(universe_ids.shape[2]):
universes[x, y, z] = \
self.get_universe_by_id(universe_ids[x, y, z])
for x in range(universe_ids.shape[2]):
universes[z, y, x] = \
self.get_universe_by_id(universe_ids[z, y, x])
# Transpose, reverse y-dimension for appropriate ordering
shape = universes.shape
universes = np.transpose(universes, (1, 0, 2))
universes.shape = shape
universes = universes[:, ::-1, :]
# Use 2D NumPy array to store lattice universes for 2D lattices
if len(dimension) == 2:
universes = np.squeeze(universes)
universes = np.atleast_2d(universes)
# Set the universes for the lattice
lattice.universes = universes
if offsets is not None:
offsets = np.swapaxes(offsets, 0, 1)
offsets = np.swapaxes(offsets, 1, 2)
offsets = np.swapaxes(offsets, 0, 2)
lattice.offsets = offsets
# Add the Lattice to the global dictionary of all Lattices
@ -474,10 +478,14 @@ class Summary(object):
# Retrieve the object corresponding to the fill type and ID
if fill_type == 'normal':
if fill_id > 0:
fill = self.get_material_by_id(fill_id)
if isinstance(fill_id, Iterable):
fill = [self.get_material_by_id(mat) if mat > 0 else 'void'
for mat in fill_id]
else:
fill = 'void'
if fill_id > 0:
fill = self.get_material_by_id(fill_id)
else:
fill = 'void'
elif fill_type == 'universe':
fill = self.get_universe_by_id(fill_id)
else:
@ -520,12 +528,18 @@ class Summary(object):
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id, tally_name)
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = self._f['{0}/moment_orders'.format(subbase)].value
# Read score metadata
scores = self._f['{0}/score_bins'.format(subbase)].value
for score in scores:
tally.add_score(score.decode())
num_score_bins = self._f['{0}/n_score_bins'.format(subbase)][...]
tally.num_score_bins = num_score_bins
for j, score in enumerate(scores):
score = score.decode()
# If this is a moment, use generic moment order
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
tally.add_score(score)
# Read filter metadata
num_filters = self._f['{0}/n_filters'.format(subbase)].value
@ -542,11 +556,11 @@ class Summary(object):
bins = self._f['{0}/bins'.format(subsubbase)][...]
# Create Filter object
filter = openmc.Filter(filter_type, bins)
filter.num_bins = num_bins
new_filter = openmc.Filter(filter_type, bins)
new_filter.num_bins = num_bins
# Add Filter to the Tally
tally.add_filter(filter)
tally.add_filter(new_filter)
# Add Tally to the global dictionary of all Tallies
self.tallies[tally_id] = tally

1243
openmc/tallies.py
View file

File diff suppressed because it is too large Load diff

58
openmc/universe.py
View file

@ -50,7 +50,7 @@ class Cell(object):
Unique identifier for the cell
name : str
Name of the cell
fill : Material or Universe or Lattice or 'void'
fill : Material or Universe or Lattice or 'void' or iterable of Material
Indicates what the region of space is filled with
region : openmc.region.Region
Region of space that is assigned to the cell.
@ -63,6 +63,8 @@ class Cell(object):
that is used to translate (shift) the universe.
offsets : ndarray
Array of offsets used for distributed cell searches
distribcell_index : int
Index of this cell in distribcell arrays
"""
@ -76,6 +78,7 @@ class Cell(object):
self._rotation = None
self._translation = None
self._offsets = None
self._distribcell_index = None
def __eq__(self, other):
if not isinstance(other, Cell):
@ -109,6 +112,12 @@ class Cell(object):
if isinstance(self._fill, openmc.Material):
string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t',
self._fill._id)
elif isinstance(self._fill, Iterable):
string += '{0: <16}{1}'.format('\tMaterial', '=\t')
string += '['
string += ', '.join(['void' if m == 'void' else str(m.id)
for m in self.fill])
string += ']\n'
elif isinstance(self._fill, (Universe, Lattice)):
string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t',
self._fill._id)
@ -122,6 +131,8 @@ class Cell(object):
string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t',
self._translation)
string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets)
string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t',
self._distribcell_index)
return string
@ -164,6 +175,10 @@ class Cell(object):
def offsets(self):
return self._offsets
@property
def distribcell_index(self):
return self._distribcell_index
@id.setter
def id(self, cell_id):
if cell_id is None:
@ -196,6 +211,11 @@ class Cell(object):
elif isinstance(fill, openmc.Material):
self._type = 'normal'
elif isinstance(fill, Iterable):
cv.check_type('cell.fill', fill, Iterable,
(openmc.Material, basestring))
self._type = 'normal'
elif isinstance(fill, Universe):
self._type = 'fill'
@ -231,6 +251,11 @@ class Cell(object):
cv.check_type('cell region', region, Region)
self._region = region
@distribcell_index.setter
def distribcell_index(self, ind):
cv.check_type('distribcell index', ind, Integral)
self._distribcell_index = ind
def add_surface(self, surface, halfspace):
"""Add a half-space to the list of half-spaces whose intersection defines the
cell.
@ -271,7 +296,7 @@ class Cell(object):
else:
self.region = Intersection(self.region, region)
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
cell_id = path[0]
path = path[1:]
@ -282,12 +307,12 @@ class Cell(object):
# If the Cell is filled by a Universe
elif self._type == 'fill':
offset = self._offsets[filter_offset-1]
offset += self._fill.get_offset(path, filter_offset)
offset = self.offsets[distribcell_index-1]
offset += self.fill.get_cell_instance(path, distribcell_index)
# If the Cell is filled by a Lattice
else:
offset = self._fill.get_offset(path, filter_offset)
offset = self.fill.get_cell_instance(path, distribcell_index)
return offset
@ -380,6 +405,10 @@ class Cell(object):
if isinstance(self._fill, openmc.Material):
element.set("material", str(self._fill._id))
elif isinstance(self._fill, Iterable):
element.set("material", ' '.join([m if m == 'void' else str(m.id)
for m in self.fill]))
elif isinstance(self._fill, (Universe, Lattice)):
element.set("fill", str(self._fill._id))
self._fill.create_xml_subelement(xml_element)
@ -591,7 +620,7 @@ class Universe(object):
self._cells.clear()
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
path = path[1:]
@ -599,7 +628,7 @@ class Universe(object):
cell_id = path[0]
# Make a recursive call to the Cell within this Universe
offset = self._cells[cell_id].get_offset(path, filter_offset)
offset = self.cells[cell_id].get_cell_instance(path, distribcell_index)
# Return the offset computed at all nested Universe levels
return offset
@ -807,7 +836,7 @@ class Lattice(object):
def universes(self, universes):
cv.check_iterable_type('lattice universes', universes, Universe,
min_depth=2, max_depth=3)
self._universes = universes
self._universes = np.asarray(universes)
def get_unique_universes(self):
"""Determine all unique universes in the lattice
@ -1059,21 +1088,22 @@ class RectLattice(Lattice):
cv.check_greater_than('lattice pitch', dim, 0.0)
self._pitch = pitch
def get_offset(self, path, filter_offset):
def get_cell_instance(self, path, distribcell_index):
# Get the current element and remove it from the list
i = path[0]
path = path[1:]
# For 2D Lattices
if len(self._dimension) == 2:
offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1]
offset += self._universes[i[1]][i[2]].get_offset(path, filter_offset)
offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1]
offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path,
distribcell_index)
# For 3D Lattices
else:
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1]
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_offset(path,
filter_offset)
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1]
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance(
path, distribcell_index)
return offset

3
setup.py
View file

@ -11,7 +11,7 @@ except ImportError:
kwargs = {'name': 'openmc',
'version': '0.7.1',
'packages': ['openmc', 'openmc.mgxs'],
'packages': ['openmc', 'openmc.mgxs', 'openmc.stats'],
'scripts': glob.glob('scripts/openmc-*'),
# Metadata
@ -37,6 +37,7 @@ if have_setuptools:
# Optional dependencies
'extras_require': {
'pandas': ['pandas'],
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],
'validate': ['lxml']
}})

4
src/ace_header.F90
View file

@ -48,7 +48,7 @@ module ace_header
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
@ -308,6 +308,8 @@ module ace_header
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)

27
src/constants.F90
View file

@ -14,7 +14,7 @@ module constants
integer, parameter :: REVISION_STATEPOINT = 14
integer, parameter :: REVISION_PARTICLE_RESTART = 1
integer, parameter :: REVISION_TRACK = 1
integer, parameter :: REVISION_SUMMARY = 1
integer, parameter :: REVISION_SUMMARY = 2
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -139,6 +139,9 @@ module constants
! Maximum number of lost particles
integer, parameter :: MAX_LOST_PARTICLES = 10
! Maximum number of lost particles, relative to the total number of particles
real(8), parameter :: REL_MAX_LOST_PARTICLES = 1e-6_8
! ============================================================================
! CROSS SECTION RELATED CONSTANTS
@ -363,28 +366,6 @@ module constants
integer, parameter :: STREAM_TALLIES = 2
integer, parameter :: STREAM_SOURCE = 3
! ============================================================================
! EXTERNAL SOURCE PARAMETERS
! Source spatial distribution types
integer, parameter :: &
SRC_SPACE_BOX = 1, & ! Source in a rectangular prism
SRC_SPACE_POINT = 2, & ! Source at a single point
SRC_SPACE_FISSION = 3 ! Source in prism filtered by fissionable mats
! Source angular distribution types
integer, parameter :: &
SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular
SRC_ANGLE_MONO = 2, & ! Monodirectional source
SRC_ANGLE_TABULAR = 3 ! Tabular distribution
! Source energy distribution types
integer, parameter :: &
SRC_ENERGY_MONO = 1, & ! Monoenergetic source
SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum
SRC_ENERGY_WATT = 3, & ! Watt fission spectrum
SRC_ENERGY_TABULAR = 4 ! Tabular distribution
! ============================================================================
! MISCELLANEOUS CONSTANTS

162
src/distribution_multivariate.F90 Normal file
View file

@ -0,0 +1,162 @@
module distribution_multivariate
use constants, only: ONE, TWO, PI
use distribution_univariate, only: Distribution
use math, only: rotate_angle
use random_lcg, only: prn
implicit none
!===============================================================================
! UNITSPHEREDISTRIBUTION type defines a probability density function for points
! on the unit sphere. Extensions of this type are used to sample angular
! distributions for starting sources
!===============================================================================
type, abstract :: UnitSphereDistribution
real(8) :: reference_uvw(3)
contains
procedure(iSample), deferred :: sample
end type UnitSphereDistribution
abstract interface
function iSample(this) result(uvw)
import UnitSphereDistribution
class(UnitSphereDistribution), intent(in) :: this
real(8) :: uvw(3)
end function iSample
end interface
!===============================================================================
! Derived classes of UnitSphereDistribution
!===============================================================================
! Explicit distribution of polar and azimuthal angles
type, extends(UnitSphereDistribution) :: PolarAzimuthal
class(Distribution), allocatable :: mu
class(Distribution), allocatable :: phi
contains
procedure :: sample => polar_azimuthal_sample
end type PolarAzimuthal
! Uniform distribution on the unit sphere
type, extends(UnitSphereDistribution) :: Isotropic
contains
procedure :: sample => isotropic_sample
end type Isotropic
! Monodirectional distribution
type, extends(UnitSphereDistribution) :: Monodirectional
contains
procedure :: sample => monodirectional_sample
end type Monodirectional
!===============================================================================
! SPATIALDISTRIBUTION type defines a probability density function for arbitrary
! points in Euclidean space.
!===============================================================================
type, abstract :: SpatialDistribution
contains
procedure(iSampleSpatial), deferred :: sample
end type SpatialDistribution
abstract interface
function iSampleSpatial(this) result(xyz)
import SpatialDistribution
class(SpatialDistribution), intent(in) :: this
real(8) :: xyz(3)
end function iSampleSpatial
end interface
type, extends(SpatialDistribution) :: CartesianIndependent
class(Distribution), allocatable :: x
class(Distribution), allocatable :: y
class(Distribution), allocatable :: z
contains
procedure :: sample => cartesian_independent_sample
end type CartesianIndependent
type, extends(SpatialDistribution) :: SpatialBox
real(8) :: lower_left(3)
real(8) :: upper_right(3)
logical :: only_fissionable = .false.
contains
procedure :: sample => spatial_box_sample
end type SpatialBox
type, extends(SpatialDistribution) :: SpatialPoint
real(8) :: xyz(3)
contains
procedure :: sample => spatial_point_sample
end type SpatialPoint
contains
function polar_azimuthal_sample(this) result(uvw)
class(PolarAzimuthal), intent(in) :: this
real(8) :: uvw(3)
real(8) :: mu ! cosine of polar angle
real(8) :: phi ! azimuthal angle
! Sample cosine of polar angle
mu = this % mu % sample()
if (mu == ONE) then
uvw(:) = this % reference_uvw
else
! Sample azimuthal angle
phi = this % phi % sample()
uvw(:) = rotate_angle(this % reference_uvw, mu, phi)
end if
end function polar_azimuthal_sample
function isotropic_sample(this) result(uvw)
class(Isotropic), intent(in) :: this
real(8) :: uvw(3)
real(8) :: phi
real(8) :: mu
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
end function isotropic_sample
function monodirectional_sample(this) result(uvw)
class(Monodirectional), intent(in) :: this
real(8) :: uvw(3)
uvw(:) = this % reference_uvw
end function monodirectional_sample
function cartesian_independent_sample(this) result(xyz)
class(CartesianIndependent), intent(in) :: this
real(8) :: xyz(3)
xyz(1) = this % x % sample()
xyz(2) = this % y % sample()
xyz(3) = this % z % sample()
end function cartesian_independent_sample
function spatial_box_sample(this) result(xyz)
class(SpatialBox), intent(in) :: this
real(8) :: xyz(3)
integer :: i
real(8) :: r(3)
r = [ (prn(), i = 1,3) ]
xyz(:) = this % lower_left + r*(this % upper_right - this % lower_left)
end function spatial_box_sample
function spatial_point_sample(this) result(xyz)
class(SpatialPoint), intent(in) :: this
real(8) :: xyz(3)
xyz(:) = this % xyz
end function spatial_point_sample
end module distribution_multivariate

347
src/distribution_univariate.F90 Normal file
View file

@ -0,0 +1,347 @@
module distribution_univariate
use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR, MAX_LINE_LEN, &
MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use string, only: to_lower
use xml_interface
implicit none
!===============================================================================
! DISTRIBUTION type defines a probability density function
!===============================================================================
type, abstract :: Distribution
contains
procedure(iSample), deferred :: sample
end type Distribution
type DistributionContainer
class(Distribution), allocatable :: obj
end type DistributionContainer
abstract interface
function iSample(this) result(x)
import Distribution
class(Distribution), intent(in) :: this
real(8) :: x
end function iSample
end interface
!===============================================================================
! Derived classes of Distribution
!===============================================================================
! Discrete distribution
type, extends(Distribution) :: Discrete
real(8), allocatable :: x(:)
real(8), allocatable :: p(:)
contains
procedure :: sample => discrete_sample
procedure :: initialize => discrete_initialize
end type Discrete
! Uniform distribution over the interval [a,b]
type, extends(Distribution) :: Uniform
real(8) :: a
real(8) :: b
contains
procedure :: sample => uniform_sample
end type Uniform
! Maxwellian distribution of form c*E*exp(-E/a)
type, extends(Distribution) :: Maxwell
real(8) :: theta
contains
procedure :: sample => maxwell_sample
end type Maxwell
! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E))
type, extends(Distribution) :: Watt
real(8) :: a
real(8) :: b
contains
procedure :: sample => watt_sample
end type Watt
! Histogram or linear-linear interpolated tabular distribution
type, extends(Distribution) :: Tabular
integer :: interpolation
real(8), allocatable :: x(:) ! tabulated independent variable
real(8), allocatable :: p(:) ! tabulated probability density
real(8), allocatable :: c(:) ! cumulative distribution at tabulated values
contains
procedure :: sample => tabular_sample
procedure :: initialize => tabular_initialize
end type Tabular
contains
function discrete_sample(this) result(x)
class(Discrete), intent(in) :: this
real(8) :: x
integer :: i ! loop counter
integer :: n ! size of distribution
real(8) :: c ! cumulative frequency
real(8) :: xi ! sampled CDF value
n = size(this%x)
if (n > 1) then
xi = prn()
c = ZERO
do i = 1, size(this%x)
c = c + this%p(i)
if (xi < c) exit
end do
x = this%x(i)
else
x = this%x(1)
end if
end function discrete_sample
subroutine discrete_initialize(this, x, p)
class(Discrete), intent(inout) :: this
real(8), intent(in) :: x(:)
real(8), intent(in) :: p(:)
integer :: n
! Check length of x, p arrays
if (size(x) /= size(p)) then
call fatal_error('Tabulated probabilities not of same length as &
&independent variable.')
end if
! Copy probability density function
n = size(x)
allocate(this%x(n), this%p(n))
this%x(:) = x(:)
this%p(:) = p(:)
! Normalize density function
this%p(:) = this%p(:)/sum(this%p)
end subroutine
function uniform_sample(this) result(x)
class(Uniform), intent(in) :: this
real(8) :: x
x = this%a + prn()*(this%b - this%a)
end function uniform_sample
function maxwell_sample(this) result(x)
class(Maxwell), intent(in) :: this
real(8) :: x
x = maxwell_spectrum(this%theta)
end function maxwell_sample
function watt_sample(this) result(x)
class(Watt), intent(in) :: this
real(8) :: x
x = watt_spectrum(this%a, this%b)
end function watt_sample
function tabular_sample(this) result(x)
class(Tabular), intent(in) :: this
real(8) :: x
integer :: i
real(8) :: c ! sampled cumulative frequency
real(8) :: m ! slope of PDF
real(8) :: x_i, x_i1 ! i-th and (i+1)th x values
real(8) :: c_i, c_i1 ! i-th and (i+1)th cumulative distribution values
real(8) :: p_i, p_i1 ! i-th and (i+1)th probability density values
! Sample value of CDF
c = prn()
! Find first CDF bin which is above the sampled value
c_i = this%c(1)
do i = 1, size(this%c) - 1
c_i1 = this%c(i + 1)
if (c <= c_i1) exit
c_i = c_i1
end do
! Determine bounding PDF values
x_i = this%x(i)
p_i = this%p(i)
if (this%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_i > ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i
end if
else
! Linear-linear interpolation
x_i1 = this%x(i + 1)
p_i1 = this%p(i + 1)
m = (p_i1 - p_i)/(x_i1 - x_i)
if (m == ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i + (sqrt(max(ZERO, p_i*p_i + 2*m*(c - c_i))) - p_i)/m
end if
end if
end function tabular_sample
subroutine tabular_initialize(this, x, p, interp)
class(Tabular), intent(inout) :: this
real(8), intent(in) :: x(:)
real(8), intent(in) :: p(:)
integer, intent(in) :: interp
integer :: i
integer :: n
! Check interpolation parameter
if (interp /= HISTOGRAM .and. interp /= LINEAR_LINEAR) then
call fatal_error('Only histogram and linear-linear interpolation for tabular &
&distribution is supported.')
end if
! Check length of x, p arrays
if (size(x) /= size(p)) then
call fatal_error('Tabulated probabilities not of same length as &
&independent variable.')
end if
! Copy probability density function and interpolation parameter
n = size(x)
allocate(this%x(n), this%p(n), this%c(n))
this%interpolation = interp
this%x(:) = x(:)
this%p(:) = p(:)
! Calculate cumulative distribution function
this%c(1) = ZERO
do i = 2, n
if (this%interpolation == HISTOGRAM) then
this%c(i) = this%c(i-1) + this%p(i-1)*(this%x(i) - this%x(i-1))
elseif (this%interpolation == LINEAR_LINEAR) then
this%c(i) = this%c(i-1) + HALF*(this%p(i-1) + this%p(i)) * &
(this%x(i) - this%x(i-1))
end if
end do
! Normalize density and distribution functions
this%p(:) = this%p(:)/this%c(n)
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
subroutine distribution_from_xml(dist, node_dist)
class(Distribution), allocatable, intent(inout) :: dist
type(Node), pointer :: node_dist
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: temp_str
integer :: n
integer :: temp_int
real(8), allocatable :: temp_real(:)
if (check_for_node(node_dist, "type")) then
! Determine type of distribution
call get_node_value(node_dist, "type", type)
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
else
n = 0
end if
! Allocate extension of Distribution
select case (to_lower(type))
case ('uniform')
allocate(Uniform :: dist)
if (n /= 2) then
call fatal_error('Uniform distribution must have two &
&parameters specified.')
end if
case ('maxwell')
allocate(Maxwell :: dist)
if (n /= 1) then
call fatal_error('Maxwell energy distribution must have one &
&parameter specified.')
end if
case ('watt')
allocate(Watt :: dist)
if (n /= 2) then
call fatal_error('Watt energy distribution must have two &
&parameters specified.')
end if
case ('discrete')
allocate(Discrete :: dist)
case ('tabular')
allocate(Tabular :: dist)
case default
call fatal_error('Invalid distribution type: ' // trim(type) // '.')
end select
! Read parameters and interpolation for distribution
select type (dist)
type is (Uniform)
allocate(temp_real(2))
call get_node_array(node_dist, "parameters", temp_real)
dist%a = temp_real(1)
dist%b = temp_real(2)
deallocate(temp_real)
type is (Maxwell)
call get_node_value(node_dist, "parameters", dist%theta)
type is (Watt)
allocate(temp_real(2))
call get_node_array(node_dist, "parameters", temp_real)
dist%a = temp_real(1)
dist%b = temp_real(2)
deallocate(temp_real)
type is (Discrete)
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n))
deallocate(temp_real)
type is (Tabular)
! Read interpolation
if (check_for_node(node_dist, "interpolation")) then
call get_node_value(node_dist, "interpolation", temp_str)
select case(to_lower(temp_str))
case ('histogram')
temp_int = HISTOGRAM
case ('linear-linear')
temp_int = LINEAR_LINEAR
case default
call fatal_error("Unknown interpolation type for distribution: " &
// trim(temp_str))
end select
else
temp_int = HISTOGRAM
end if
! Read and initialize tabular distribution
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n), temp_int)
deallocate(temp_real)
end select
end if
end subroutine distribution_from_xml
end module distribution_univariate

15
src/error.F90
View file

@ -3,8 +3,6 @@ module error
use, intrinsic :: ISO_FORTRAN_ENV
use constants
use global
#ifdef MPI
use message_passing
#endif
@ -87,6 +85,9 @@ contains
integer :: line_wrap ! length of line
integer :: length ! length of message
integer :: indent ! length of indentation
#ifdef MPI
integer :: mpi_err
#endif
! set default error code
@ -136,16 +137,6 @@ contains
end if
end do
! Write information on current batch, generation, and particle
if (current_batch > 0) then
write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch
write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen
write(ERROR_UNIT,*)
end if
! Release memory from all allocatable arrays
call free_memory()
#ifdef MPI
! Abort MPI
call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err)

58
src/geometry.F90
View file

@ -193,7 +193,9 @@ contains
logical, intent(inout) :: found
integer, optional :: search_cells(:)
integer :: i ! index over cells
integer :: j ! coordinate level index
integer :: j, k ! coordinate level index
integer :: offset ! instance # of a distributed cell
integer :: distribcell_index
integer :: i_xyz(3) ! indices in lattice
integer :: n ! number of cells to search
integer :: index_cell ! index in cells array
@ -246,9 +248,36 @@ contains
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
! set material
! Set the particle material
p % last_material = p % material
p % material = c % material
if (size(c % material) == 1) then
! Only one material for this cell; assign that one to the particle.
p % material = c % material(1)
else
! Distributed instances of this cell have different materials.
! Determine which instance this is and assign the matching material.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
end if
end if
end do
p % material = c % material(offset + 1)
end if
elseif (c % type == CELL_FILL) then CELL_TYPE
! ======================================================================
@ -308,9 +337,10 @@ contains
else
! Particle is outside the lattice.
if (lat % outer == NO_OUTER_UNIVERSE) then
call fatal_error("A particle is outside latttice " &
// trim(to_str(lat % id)) // " but the lattice has no &
&defined outer universe.")
call handle_lost_particle(p, "Particle " // trim(to_str(p %id)) &
// " is outside lattice " // trim(to_str(lat % id)) &
// " but the lattice has no defined outer universe.")
return
else
p % coord(j + 1) % universe = lat % outer
end if
@ -536,9 +566,11 @@ contains
p % n_coord = 1
call find_cell(p, found)
if (.not. found) then
call handle_lost_particle(p, "Could not locate particle " &
// trim(to_str(p % id)) // " after crossing a lattice boundary.")
return
if (p % alive) then ! Particle may have been killed in find_cell
call handle_lost_particle(p, "Could not locate particle " &
// trim(to_str(p % id)) // " after crossing a lattice boundary.")
return
end if
end if
else OUTSIDE_LAT
@ -947,6 +979,8 @@ contains
type(Particle), intent(inout) :: p
character(*) :: message
integer(8) :: tot_n_particles
! Print warning and write lost particle file
call warning(message)
call write_particle_restart(p)
@ -956,9 +990,13 @@ contains
!$omp atomic
n_lost_particles = n_lost_particles + 1
! Count the total number of simulated particles (on this processor)
tot_n_particles = n_batches * gen_per_batch * work
! Abort the simulation if the maximum number of lost particles has been
! reached
if (n_lost_particles == MAX_LOST_PARTICLES) then
if (n_lost_particles >= MAX_LOST_PARTICLES .and. &
n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then
call fatal_error("Maximum number of lost particles has been reached.")
end if

16
src/geometry_header.F90
View file

@ -126,15 +126,19 @@ module geometry_header
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in
! the geom
integer :: material ! Material within cell (0 for
! universe)
integer, allocatable :: offset (:) ! Distribcell offset for tally
integer, allocatable :: material(:) ! Material within cell. Multiple
! materials for distribcell
! instances. 0 signifies a universe
integer, allocatable :: offset(:) ! Distribcell offset for tally
! counter
integer, allocatable :: region(:) ! Definition of spatial region as
! Boolean expression of half-spaces
! Boolean expression of half-spaces
integer, allocatable :: rpn(:) ! Reverse Polish notation for region
! expression
logical :: simple ! Is the region simple (intersections only)
! expression
logical :: simple ! Is the region simple (intersections
! only)
integer :: distribcell_index ! Index corresponding to this cell in
! distribcell arrays
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

14
src/global.F90
View file

@ -12,7 +12,7 @@ module global
use plot_header, only: ObjectPlot
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: ExtSource
use source_header, only: SourceDistribution
use tally_header, only: TallyObject, TallyMap, TallyResult
use trigger_header, only: KTrigger
use timer_header, only: Timer
@ -182,7 +182,7 @@ module global
logical :: satisfy_triggers = .false. ! whether triggers are satisfied
! External source
type(ExtSource), target :: external_source
type(SourceDistribution), allocatable :: external_source(:)
! Source and fission bank
type(Bank), allocatable, target :: source_bank(:)
@ -289,9 +289,6 @@ module global
character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart
character(MAX_FILE_LEN) :: path_output = '' ! Path to output directory
! Random number seed
integer(8) :: seed = 1_8
! The verbosity controls how much information will be printed to the
! screen and in logs
integer :: verbosity = 7
@ -453,12 +450,7 @@ contains
if (allocated(micro_xs)) deallocate(micro_xs)
! Deallocate external source
if (allocated(external_source % params_space)) &
deallocate(external_source % params_space)
if (allocated(external_source % params_angle)) &
deallocate(external_source % params_angle)
if (allocated(external_source % params_energy)) &
deallocate(external_source % params_energy)
if (allocated(external_source)) deallocate(external_source)
! Deallocate k and entropy
if (allocated(k_generation)) deallocate(k_generation)

2
src/hdf5_interface.F90
View file

@ -1898,7 +1898,7 @@ contains
logical :: mpio
integer :: hdf5_err
integer :: driver
integer(HID_T) :: driver
integer(HID_T) :: file_id
integer(HID_T) :: fapl_id

254
src/initialize.F90
View file

@ -614,31 +614,33 @@ contains
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
id = c%material
if (id == MATERIAL_VOID) then
c%type = CELL_NORMAL
elseif (id /= 0) then
if (material_dict%has_key(id)) then
c%type = CELL_NORMAL
c%material = material_dict%get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
&// " specified on cell " // trim(to_str(c%id)))
end if
else
id = c%fill
if (universe_dict%has_key(id)) then
c%type = CELL_FILL
c%fill = universe_dict%get_key(id)
elseif (lattice_dict%has_key(id)) then
lid = lattice_dict%get_key(id)
c%type = CELL_LATTICE
c%fill = lid
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has_key(id)) then
c % type = CELL_FILL
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % type = CELL_LATTICE
c % fill = lid
else
call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "&
&// trim(to_str(c%id)) // " is neither a universe nor a &
// trim(to_str(c % id)) // " is neither a universe nor a &
&lattice.")
end if
else
do j = 1, size(c % material)
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = CELL_NORMAL
else if (material_dict % has_key(id)) then
c % type = CELL_NORMAL
c % material(j) = material_dict % get_key(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
// " specified on cell " // trim(to_str(c % id)))
end if
end do
end if
end do
@ -952,86 +954,82 @@ contains
subroutine prepare_distribcell()
integer :: i, j ! Tally, filter loop counters
integer :: n_filt ! Number of filters originally in tally
logical :: count_all ! Count all cells
type(TallyObject), pointer :: t ! Current tally
type(Universe), pointer :: univ ! Pointer to universe
type(Cell), pointer :: c ! Pointer to cell
integer :: i, j ! Tally, filter loop counters
logical :: distribcell_active ! Does simulation use distribcell?
integer, allocatable :: univ_list(:) ! Target offsets
integer, allocatable :: counts(:,:) ! Target count
logical, allocatable :: found(:,:) ! Target found
count_all = .false.
! Assume distribcell is not needed until proven otherwise.
distribcell_active = .false.
! Loop over tallies
! We need distribcell if any tallies have distribcell filters.
do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
n_filt = t%n_filters
! Loop over the filters to determine how many additional filters
! need to be added to this tally
do j = 1, t%n_filters
! Determine type of filter
if (t%filters(j)%type == FILTER_DISTRIBCELL) then
count_all = .true.
if (size(t%filters(j)%int_bins) > 1) then
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
distribcell_active = .true.
if (size(tallies(i) % filters(j) % int_bins) > 1) then
call fatal_error("A distribcell filter was specified with &
&multiple bins. This feature is not supported.")
end if
end if
end do
end do
if (count_all) then
univ => universes(BASE_UNIVERSE)
! sum the number of occurrences of all cells
call count_instance(univ)
! Loop over tallies
do i = 1, n_tallies
! Get pointer to tally
t => tallies(i)
! Initialize the filters
do j = 1, t%n_filters
! Set the number of bins to the number of instances of the cell
if (t%filters(j)%type == FILTER_DISTRIBCELL) then
c => cells(t%filters(j)%int_bins(1))
t%filters(j)%n_bins = c%instances
end if
end do
! We also need distribcell if any distributed materials are present.
if (.not. distribcell_active) then
do i = 1, n_cells
if (size(cells(i) % material) > 1) then
distribcell_active = .true.
exit
end if
end do
end if
! If distribcell isn't used in this simulation then no more work left to do.
if (.not. distribcell_active) return
! Count the number of instances of each cell.
call count_instance(universes(BASE_UNIVERSE))
! Set the number of bins in all distribcell filters.
do i = 1, n_tallies
do j = 1, tallies(i) % n_filters
associate (filt => tallies(i) % filters(j))
if (filt % type == FILTER_DISTRIBCELL) then
! Set the number of bins to the number of instances of the cell.
filt % n_bins = cells(filt % int_bins(1)) % instances
end if
end associate
end do
end do
! Make sure the number of materials matches the number of cell instances for
! distributed materials.
do i = 1, n_cells
associate (c => cells(i))
if (size(c % material) > 1) then
if (size(c % material) /= c % instances) then
call fatal_error("Cell " // trim(to_str(c % id)) // " was &
&specified with " // trim(to_str(size(c % material))) &
// " materials but has " // trim(to_str(c % instances)) &
// " distributed instances. The number of materials must &
&equal one or the number of instances.")
end if
end if
end associate
end do
! Allocate offset maps at each level in the geometry
call allocate_offsets(univ_list, counts, found)
! Calculate offsets for each target distribcell
do i = 1, n_maps
do j = 1, n_universes
univ => universes(j)
call calc_offsets(univ_list(i), i, univ, counts, found)
call calc_offsets(univ_list(i), i, universes(j), counts, found)
end do
end do
! Deallocate temporary target variable arrays
deallocate(counts)
deallocate(found)
deallocate(univ_list)
end subroutine prepare_distribcell
!===============================================================================
@ -1045,38 +1043,33 @@ contains
integer, intent(out), allocatable :: counts(:,:) ! Target count
logical, intent(out), allocatable :: found(:,:) ! Target found
integer :: i, j, k, l, m ! Loop counters
type(SetInt) :: cell_list ! distribells to track
type(Universe), pointer :: univ ! pointer to universe
class(Lattice), pointer :: lat ! pointer to lattice
type(TallyObject), pointer :: t ! pointer to tally
type(TallyFilter), pointer :: filter ! pointer to filter
integer :: i, j, k ! Loop counters
type(SetInt) :: cell_list ! distribells to track
! Begin gathering list of cells in distribcell tallies
n_maps = 0
! Populate list of distribcells to track
! List all cells referenced in distribcell filters.
do i = 1, n_tallies
t => tallies(i)
do j = 1, t%n_filters
filter => t%filters(j)
if (filter%type == FILTER_DISTRIBCELL) then
if (.not. cell_list%contains(filter%int_bins(1))) then
call cell_list%add(filter%int_bins(1))
end if
do j = 1, tallies(i) % n_filters
if (tallies(i) % filters(j) % type == FILTER_DISTRIBCELL) then
call cell_list % add(tallies(i) % filters(j) % int_bins(1))
end if
end do
end do
! List all cells with multiple (distributed) materials.
do i = 1, n_cells
if (size(cells(i) % material) > 1) then
call cell_list % add(i)
end if
end do
! Compute the number of unique universes containing these distribcells
! to determine the number of offset tables to allocate
do i = 1, n_universes
univ => universes(i)
do j = 1, univ%n_cells
if (cell_list%contains(univ%cells(j))) then
do j = 1, universes(i) % n_cells
if (cell_list % contains(universes(i) % cells(j))) then
n_maps = n_maps + 1
end if
end do
@ -1085,42 +1078,23 @@ contains
! Allocate the list of offset tables for each unique universe
allocate(univ_list(n_maps))
! Allocate list to accumulate target distribccell counts in each universe
! Allocate list to accumulate target distribcell counts in each universe
allocate(counts(n_universes, n_maps))
counts(:,:) = 0
! Allocate list to track if target distribcells are found in each universe
allocate(found(n_universes, n_maps))
counts(:,:) = 0
found(:,:) = .false.
! Search through universes for distributed cells and assign each one a
! unique distribcell array index.
k = 1
do i = 1, n_universes
univ => universes(i)
do j = 1, univ%n_cells
if (cell_list%contains(univ%cells(j))) then
! Loop over all tallies
do l = 1, n_tallies
t => tallies(l)
do m = 1, t%n_filters
filter => t%filters(m)
! Loop over only distribcell filters
! If filter points to cell we just found, set offset index
if (filter%type == FILTER_DISTRIBCELL) then
if (filter%int_bins(1) == univ%cells(j)) then
filter%offset = k
end if
end if
end do
end do
univ_list(k) = univ%id
do j = 1, universes(i) % n_cells
if (cell_list % contains(universes(i) % cells(j))) then
cells(universes(i) % cells(j)) % distribcell_index = k
univ_list(k) = universes(i) % id
k = k + 1
end if
end do
@ -1128,29 +1102,31 @@ contains
! Allocate the offset tables for lattices
do i = 1, n_lattices
lat => lattices(i)%obj
associate(lat => lattices(i) % obj)
select type(lat)
select type(lat)
type is (RectLattice)
allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), &
lat%n_cells(3)))
type is (HexLattice)
allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, &
2 * lat%n_rings - 1, lat%n_axial))
end select
lat%offset(:, :, :, :) = 0
type is (RectLattice)
allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), &
lat % n_cells(3)))
type is (HexLattice)
allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, &
2 * lat % n_rings - 1, lat % n_axial))
end select
lat % offset(:, :, :, :) = 0
end associate
end do
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i)%material == NONE) then
allocate(cells(i)%offset(n_maps))
if (cells(i) % type /= CELL_NORMAL) then
allocate(cells(i) % offset(n_maps))
end if
end do
! Free up memory
call cell_list % clear()
end subroutine allocate_offsets
end module initialize

448
src/input_xml.F90
View file

@ -3,6 +3,8 @@ module input_xml
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
use distribution_multivariate
use distribution_univariate
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
@ -11,7 +13,7 @@ module input_xml
use mesh_header, only: RegularMesh
use output, only: write_message
use plot_header
use random_lcg, only: prn
use random_lcg, only: prn, seed
use surface_header
use stl_vector, only: VectorInt
use string, only: to_lower, to_str, str_to_int, str_to_real, &
@ -54,11 +56,11 @@ contains
character(MAX_LINE_LEN) :: temp_str
integer :: i
integer :: n
integer :: coeffs_reqd
integer :: temp_int
integer :: temp_int_array3(3)
integer, allocatable :: temp_int_array(:)
integer(8) :: temp_long
real(8), allocatable :: temp_real(:)
integer :: n_tracks
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable
@ -67,6 +69,8 @@ contains
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mode => null()
type(Node), pointer :: node_source => null()
type(Node), pointer :: node_space => null()
type(Node), pointer :: node_angle => null()
type(Node), pointer :: node_dist => null()
type(Node), pointer :: node_cutoff => null()
type(Node), pointer :: node_entropy => null()
@ -79,6 +83,7 @@ contains
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_keff_trigger => null()
type(NodeList), pointer :: node_scat_list => null()
type(NodeList), pointer :: node_source_list => null()
! Display output message
call write_message("Reading settings XML file...", 5)
@ -339,185 +344,249 @@ contains
! ==========================================================================
! EXTERNAL SOURCE
! Get pointer to source
if (check_for_node(doc, "source")) then
call get_node_ptr(doc, "source", node_source)
else
call fatal_error("No source specified in settings XML file.")
end if
! Get point to list of <source> elements and make sure there is at least one
call get_node_list(doc, "source", node_source_list)
n = get_list_size(node_source_list)
if (n == 0) call fatal_error("No source specified in settings XML file.")
! Check if we want to write out source
if (check_for_node(node_source, "write_initial")) then
call get_node_value(node_source, "write_initial", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
write_initial_source = .true.
end if
! Allocate array for sources
allocate(external_source(n))
! Check for external source file
if (check_for_node(node_source, "file")) then
! Copy path of source file
call get_node_value(node_source, "file", path_source)
! Read each source
do i = 1, n
! Get pointer to source
call get_list_item(node_source_list, i, node_source)
! Check if source file exists
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Binary source file '" // trim(path_source) &
&// "' does not exist!")
! Check if we want to write out source
if (check_for_node(node_source, "write_initial")) then
call get_node_value(node_source, "write_initial", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
write_initial_source = .true.
end if
else
! Spatial distribution for external source
if (check_for_node(node_source, "space")) then
! Get pointer to spatial distribution
call get_node_ptr(node_source, "space", node_dist)
! Check for type of spatial distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('box')
external_source % type_space = SRC_SPACE_BOX
coeffs_reqd = 6
case ('fission')
external_source % type_space = SRC_SPACE_FISSION
coeffs_reqd = 6
case ('point')
external_source % type_space = SRC_SPACE_POINT
coeffs_reqd = 3
case default
call fatal_error("Invalid spatial distribution for external source: "&
&// trim(type))
end select
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
else
n = 0
end if
! Read parameters for spatial distribution
if (n < coeffs_reqd) then
call fatal_error("Not enough parameters specified for spatial &
&distribution of external source.")
elseif (n > coeffs_reqd) then
call fatal_error("Too many parameters specified for spatial &
&distribution of external source.")
elseif (n > 0) then
allocate(external_source % params_space(n))
call get_node_array(node_dist, "parameters", &
external_source % params_space)
end if
! Check for source strength
if (check_for_node(node_source, "strength")) then
call get_node_value(node_source, "strength", external_source(i)%strength)
else
call fatal_error("No spatial distribution specified for external &
&source.")
external_source(i)%strength = ONE
end if
! Determine external source angular distribution
if (check_for_node(node_source, "angle")) then
! Check for external source file
if (check_for_node(node_source, "file")) then
! Copy path of source file
call get_node_value(node_source, "file", path_source)
! Get pointer to angular distribution
call get_node_ptr(node_source, "angle", node_dist)
! Check for type of angular distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('isotropic')
external_source % type_angle = SRC_ANGLE_ISOTROPIC
coeffs_reqd = 0
case ('monodirectional')
external_source % type_angle = SRC_ANGLE_MONO
coeffs_reqd = 3
case ('tabular')
external_source % type_angle = SRC_ANGLE_TABULAR
case default
call fatal_error("Invalid angular distribution for external source: "&
&// trim(type))
end select
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
else
n = 0
! Check if source file exists
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Binary source file '" // trim(path_source) &
&// "' does not exist!")
end if
! Read parameters for angle distribution
if (n < coeffs_reqd) then
call fatal_error("Not enough parameters specified for angle &
&distribution of external source.")
elseif (n > coeffs_reqd) then
call fatal_error("Too many parameters specified for angle &
&distribution of external source.")
elseif (n > 0) then
allocate(external_source % params_angle(n))
call get_node_array(node_dist, "parameters", &
external_source % params_angle)
end if
else
! Set default angular distribution isotropic
external_source % type_angle = SRC_ANGLE_ISOTROPIC
end if
! Determine external source energy distribution
if (check_for_node(node_source, "energy")) then
! Spatial distribution for external source
if (check_for_node(node_source, "space")) then
! Get pointer to energy distribution
call get_node_ptr(node_source, "energy", node_dist)
! Get pointer to spatial distribution
call get_node_ptr(node_source, "space", node_space)
! Check for type of energy distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('monoenergetic')
external_source % type_energy = SRC_ENERGY_MONO
coeffs_reqd = 1
case ('maxwell')
external_source % type_energy = SRC_ENERGY_MAXWELL
coeffs_reqd = 1
case ('watt')
external_source % type_energy = SRC_ENERGY_WATT
coeffs_reqd = 2
case ('tabular')
external_source % type_energy = SRC_ENERGY_TABULAR
case default
call fatal_error("Invalid energy distribution for external source: " &
&// trim(type))
end select
! Check for type of spatial distribution
type = ''
if (check_for_node(node_space, "type")) &
call get_node_value(node_space, "type", type)
select case (to_lower(type))
case ('cartesian')
allocate(CartesianIndependent :: external_source(i)%space)
case ('box')
allocate(SpatialBox :: external_source(i)%space)
case ('fission')
allocate(SpatialBox :: external_source(i)%space)
select type(space => external_source(i)%space)
type is (SpatialBox)
space%only_fissionable = .true.
end select
case ('point')
allocate(SpatialPoint :: external_source(i)%space)
case default
call fatal_error("Invalid spatial distribution for external source: "&
// trim(type))
end select
select type (space => external_source(i)%space)
type is (CartesianIndependent)
! Read distribution for x coordinate
if (check_for_node(node_space, "x")) then
call get_node_ptr(node_space, "x", node_dist)
call distribution_from_xml(space%x, node_dist)
else
allocate(Discrete :: space%x)
select type (dist => space%x)
type is (Discrete)
allocate(dist%x(1), dist%p(1))
dist%x(1) = ZERO
dist%p(1) = ONE
end select
end if
! Read distribution for y coordinate
if (check_for_node(node_space, "y")) then
call get_node_ptr(node_space, "y", node_dist)
call distribution_from_xml(space%y, node_dist)
else
allocate(Discrete :: space%y)
select type (dist => space%y)
type is (Discrete)
allocate(dist%x(1), dist%p(1))
dist%x(1) = ZERO
dist%p(1) = ONE
end select
end if
if (check_for_node(node_space, "z")) then
call get_node_ptr(node_space, "z", node_dist)
call distribution_from_xml(space%z, node_dist)
else
allocate(Discrete :: space%z)
select type (dist => space%z)
type is (Discrete)
allocate(dist%x(1), dist%p(1))
dist%x(1) = ZERO
dist%p(1) = ONE
end select
end if
type is (SpatialBox)
! Make sure correct number of parameters are given
if (get_arraysize_double(node_space, "parameters") /= 6) then
call fatal_error('Box/fission spatial source must have &
&six parameters specified.')
end if
! Read lower-right/upper-left coordinates
allocate(temp_real(6))
call get_node_array(node_space, "parameters", temp_real)
space%lower_left(:) = temp_real(1:3)
space%upper_right(:) = temp_real(4:6)
deallocate(temp_real)
type is (SpatialPoint)
! Make sure correct number of parameters are given
if (get_arraysize_double(node_space, "parameters") /= 3) then
call fatal_error('Point spatial source must have &
&three parameters specified.')
end if
! Read location of point source
allocate(temp_real(3))
call get_node_array(node_space, "parameters", temp_real)
space%xyz(:) = temp_real
deallocate(temp_real)
end select
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
else
n = 0
call fatal_error("No spatial distribution specified for external &
&source.")
end if
! Read parameters for energy distribution
if (n < coeffs_reqd) then
call fatal_error("Not enough parameters specified for energy &
&distribution of external source.")
elseif (n > coeffs_reqd) then
call fatal_error("Too many parameters specified for energy &
&distribution of external source.")
elseif (n > 0) then
allocate(external_source % params_energy(n))
call get_node_array(node_dist, "parameters", &
external_source % params_energy)
! Determine external source angular distribution
if (check_for_node(node_source, "angle")) then
! Get pointer to angular distribution
call get_node_ptr(node_source, "angle", node_angle)
! Check for type of angular distribution
type = ''
if (check_for_node(node_angle, "type")) &
call get_node_value(node_angle, "type", type)
select case (to_lower(type))
case ('isotropic')
allocate(Isotropic :: external_source(i)%angle)
case ('monodirectional')
allocate(Monodirectional :: external_source(i)%angle)
case ('mu-phi')
allocate(PolarAzimuthal :: external_source(i)%angle)
case default
call fatal_error("Invalid angular distribution for external source: "&
// trim(type))
end select
! Read reference directional unit vector
if (check_for_node(node_angle, "reference_uvw")) then
n = get_arraysize_double(node_angle, "reference_uvw")
if (n /= 3) then
call fatal_error('Angular distribution reference direction must have &
&three parameters specified.')
end if
call get_node_array(node_angle, "reference_uvw", &
external_source(i)%angle%reference_uvw)
else
! By default, set reference unit vector to be positive z-direction
external_source(i)%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Read parameters for angle distribution
select type (angle => external_source(i)%angle)
type is (Monodirectional)
call get_node_array(node_angle, "reference_uvw", &
external_source(i)%angle%reference_uvw)
type is (PolarAzimuthal)
if (check_for_node(node_angle, "mu")) then
call get_node_ptr(node_angle, "mu", node_dist)
call distribution_from_xml(angle%mu, node_dist)
else
allocate(Uniform :: angle%mu)
select type (mu => angle%mu)
type is (Uniform)
mu%a = -ONE
mu%b = ONE
end select
end if
if (check_for_node(node_angle, "phi")) then
call get_node_ptr(node_angle, "phi", node_dist)
call distribution_from_xml(angle%phi, node_dist)
else
allocate(Uniform :: angle%phi)
select type (phi => angle%phi)
type is (Uniform)
phi%a = ZERO
phi%b = TWO*PI
end select
end if
end select
else
! Set default angular distribution isotropic
allocate(Isotropic :: external_source(i)%angle)
external_source(i)%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Determine external source energy distribution
if (check_for_node(node_source, "energy")) then
call get_node_ptr(node_source, "energy", node_dist)
call distribution_from_xml(external_source(i)%energy, node_dist)
else
! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
allocate(Watt :: external_source(i)%energy)
select type(energy => external_source(i)%energy)
type is (Watt)
energy%a = 0.988_8
energy%b = 2.249_8
end select
end if
else
! Set default energy distribution to Watt fission spectrum
external_source % type_energy = SRC_ENERGY_WATT
allocate(external_source % params_energy(2))
external_source % params_energy = (/ 0.988_8, 2.249_8 /)
end if
end if
end do
! Survival biasing
if (check_for_node(doc, "survival_biasing")) then
@ -982,8 +1051,7 @@ contains
subroutine read_geometry_xml()
integer :: i, j, k, m, i_x, i_a, input_index
integer :: n
integer :: n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats
integer :: universe_num
integer :: n_cells_in_univ
integer :: coeffs_reqd
@ -994,6 +1062,7 @@ contains
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
character(MAX_WORD_LEN), allocatable :: sarray(:)
character(REGION_SPEC_LEN) :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
@ -1049,8 +1118,9 @@ contains
do i = 1, n_cells
c => cells(i)
! Initialize the number of cell instances - this is a base case for distribcells
! Initialize distribcell instances and distribcell index
c % instances = 0
c % distribcell_index = NONE
! Get pointer to i-th cell node
call get_list_item(node_cell_list, i, node_cell)
@ -1085,36 +1155,50 @@ contains
end if
! Read material
word = ''
if (check_for_node(node_cell, "material")) &
call get_node_value(node_cell, "material", word)
select case(to_lower(word))
case ('void')
c % material = MATERIAL_VOID
if (check_for_node(node_cell, "material")) then
n_mats = get_arraysize_string(node_cell, "material")
case ('')
! This case is called if no material was specified
c % material = NONE
if (n_mats > 0) then
allocate(sarray(n_mats))
call get_node_array(node_cell, "material", sarray)
case default
c % material = int(str_to_int(word), 4)
allocate(c % material(n_mats))
do j = 1, n_mats
select case(trim(to_lower(sarray(j))))
case ('void')
c % material(j) = MATERIAL_VOID
case default
c % material(j) = int(str_to_int(sarray(j)), 4)
! Check for error
if (c % material == ERROR_INT) then
call fatal_error("Invalid material specified on cell " &
&// to_str(c % id))
! Check for error
if (c % material(j) == ERROR_INT) then
call fatal_error("Invalid material specified on cell " &
// to_str(c % id))
end if
end select
end do
deallocate(sarray)
else
allocate(c % material(1))
c % material(1) = NONE
end if
end select
else
allocate(c % material(1))
c % material(1) = NONE
end if
! Check to make sure that either material or fill was specified
if (c % material == NONE .and. c % fill == NONE) then
if (c % material(1) == NONE .and. c % fill == NONE) then
call fatal_error("Neither material nor fill was specified for cell " &
&// trim(to_str(c % id)))
// trim(to_str(c % id)))
end if
! Check to make sure that both material and fill haven't been
! specified simultaneously
if (c % material /= NONE .and. c % fill /= NONE) then
if (c % material(1) /= NONE .and. c % fill /= NONE) then
call fatal_error("Cannot specify material and fill simultaneously")
end if

54
src/math.F90
View file

@ -557,6 +557,60 @@ contains
end function calc_rn
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu, phi) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8), optional :: phi ! azimuthal angle
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi_ ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi) if none provided
if (present(phi)) then
phi_ = phi
else
phi_ = TWO * PI * prn()
end if
! Precompute factors to save flops
sinphi = sin(phi_)
cosphi = cos(phi_)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
!===============================================================================
! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based
! on a direct sampling scheme. The probability distribution function for a

37
src/output.F90
View file

@ -1396,8 +1396,7 @@ contains
label = ''
univ => universes(BASE_UNIVERSE)
offset = 0
call find_offset(t % filters(i_filter) % offset, &
t % filters(i_filter) % int_bins(1), &
call find_offset(t % filters(i_filter) % int_bins(1), &
univ, bin-1, offset, label)
case (FILTER_SURFACE)
i = t % filters(i_filter) % int_bins(bin)
@ -1431,15 +1430,15 @@ contains
! with the given offset
!===============================================================================
recursive subroutine find_offset(map, goal, univ, final, offset, path)
recursive subroutine find_offset(goal, univ, final, offset, path)
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
integer, intent(in) :: goal ! The target cell index
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: map ! Index in maps vector
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
@ -1454,6 +1453,9 @@ contains
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice
! Get the distribcell index for this cell
map = cells(goal) % distribcell_index
n = univ % n_cells
! Write to the geometry stack
@ -1465,17 +1467,13 @@ contains
! Look through all cells in this universe
do i = 1, n
cell_index = univ % cells(i)
c => cells(cell_index)
! If the cell ID matches the goal and the offset matches final,
! write to the geometry stack
if (cell_dict % get_key(c % id) == goal .AND. offset == final) then
path = trim(path) // "->" // to_str(c%id)
! If the cell matches the goal and the offset matches final, write to the
! geometry stack
if (univ % cells(i) == goal .AND. offset == final) then
c => cells(univ % cells(i))
path = trim(path) // "->" // to_str(c % id)
return
end if
end do
! Find the fill cell or lattice cell that we need to enter
@ -1555,7 +1553,7 @@ contains
offset = c % offset(map) + offset
next_univ => universes(c % fill)
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
! ====================================================================
@ -1595,7 +1593,7 @@ contains
path = trim(path) // "(" // trim(to_str(k)) // &
"," // trim(to_str(l)) // "," // &
trim(to_str(m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
@ -1611,7 +1609,7 @@ contains
path = trim(path) // "(" // trim(to_str(old_k)) // &
"," // trim(to_str(old_l)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
end if
@ -1656,8 +1654,7 @@ contains
trim(to_str(k - lat % n_rings)) // "," // &
trim(to_str(l - lat % n_rings)) // "," // &
trim(to_str(m)) // ")"
call find_offset(map, goal, next_univ, final, offset, &
path)
call find_offset(goal, next_univ, final, offset, path)
return
else
old_m = m
@ -1674,7 +1671,7 @@ contains
trim(to_str(old_k - lat % n_rings)) // "," // &
trim(to_str(old_l - lat % n_rings)) // "," // &
trim(to_str(old_m)) // ")"
call find_offset(map, goal, next_univ, final, offset, path)
call find_offset(goal, next_univ, final, offset, path)
return
end if

51
src/physics.F90
View file

@ -9,7 +9,7 @@ module physics
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use math, only: maxwell_spectrum, watt_spectrum
use math, only: rotate_angle, maxwell_spectrum, watt_spectrum
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: Particle
@ -1517,55 +1517,6 @@ contains
end function sample_angle
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi)
phi = TWO * PI * prn()
! Precompute factors to save flops
sinphi = sin(phi)
cosphi = cos(phi)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
!===============================================================================
! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary
! neutron from a collision or for a prompt/delayed fission neutron

14
src/plot.F90
View file

@ -82,17 +82,17 @@ contains
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord(j) % cell)
if (c % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else if (c % type == CELL_FILL) then
if (c % type == CELL_FILL) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1
else if (p % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else
rgb = pl % colors(c % material) % rgb
id = materials(c % material) % id
rgb = pl % colors(p % material) % rgb
id = materials(p % material) % id
end if
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell

5
src/random_lcg.F90
View file

@ -7,6 +7,9 @@ module random_lcg
private
save
! Random number seed
integer(8), public :: seed = 1_8
integer(8) :: prn_seed0 ! original seed
integer(8) :: prn_seed(N_STREAMS) ! current seed
integer(8) :: prn_mult ! multiplication factor, g
@ -56,8 +59,6 @@ contains
subroutine initialize_prng()
use global, only: seed
integer :: i
prn_seed0 = seed

4
src/relaxng/geometry.rnc
View file

@ -6,8 +6,8 @@ element geometry {
(element universe { xsd:int } | attribute universe { xsd:int })? &
(
(element fill { xsd:int } | attribute fill { xsd:int }) |
(element material { ( xsd:int | "void" ) } |
attribute material { ( xsd:int | "void" ) })
(element material { ( xsd:int | "void" )+ } |
attribute material { ( xsd:int | "void" )+ })
) &
(element region { xsd:string } | attribute region { xsd:string })? &
(element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &

20
src/relaxng/geometry.rng
View file

@ -47,16 +47,20 @@
</choice>
<choice>
<element name="material">
<choice>
<data type="int"/>
<value>void</value>
</choice>
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
</element>
<attribute name="material">
<choice>
<data type="int"/>
<value>void</value>
</choice>
<oneOrMore>
<choice>
<data type="int"/>
<value>void</value>
</choice>
</oneOrMore>
</attribute>
</choice>
</choice>

60
src/relaxng/settings.rnc
View file

@ -65,36 +65,36 @@ element settings {
element seed { xsd:positiveInteger }? &
element source {
element file { xsd:string { maxLength = "255" } }? &
element space {
(element type { xsd:string { maxLength = "16" } } |
attribute type { xsd:string { maxLength = "16" } }) &
(element length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { xsd:string { maxLength = "10" } })? &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })?
}? &
element angle {
(element type { xsd:string { maxLength = "16" } } |
attribute type { xsd:string { maxLength = "16" } }) &
(element length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { xsd:string { maxLength = "10" } })? &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })?
}? &
element energy {
(element type { xsd:string { maxLength = "16" } } |
attribute type { xsd:string { maxLength = "16" } }) &
(element length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { xsd:string { maxLength = "10" } })? &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })?
}? &
(element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })?
}? &
grammar {
start =
(element strength { xsd:double } | attribute strength { xsd:double })? &
(element file { xsd:string } | attribute file { xsd:string })? &
element space {
(element type { xsd:string } | attribute type { xsd:string }) &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })? &
element x { distribution }? &
element y { distribution }? &
element z { distribution }?
}? &
element angle {
(element type { xsd:string } | attribute type { xsd:string }) &
(element reference_uvw { list { xsd:double, xsd:double, xsd:double } } |
attribute reference_uvw { list { xsd:double, xsd:double, xsd:double } })? &
element mu { distribution }? &
element phi { distribution }?
}? &
element energy { distribution }? &
(element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })?
distribution =
(element type { xsd:string { maxLength = "16" } } |
attribute type { xsd:string { maxLength = "16" } }) &
(element interpolation { xsd:string } |
attribute interpolation { xsd:string })? &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })?
}
}* &
element state_point {
(

353
src/relaxng/settings.rng
View file

@ -264,209 +264,184 @@
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<zeroOrMore>
<element name="source">
<interleave>
<optional>
<element name="file">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
<optional>
<element name="space">
<interleave>
<grammar>
<start>
<interleave>
<optional>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<element name="strength">
<data type="double"/>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<attribute name="strength">
<data type="double"/>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="angle">
<interleave>
</optional>
<optional>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<element name="file">
<data type="string"/>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<attribute name="file">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="energy">
<interleave>
</optional>
<optional>
<element name="space">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<element name="x">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="y">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="z">
<ref name="distribution"/>
</element>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="angle">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="reference_uvw">
<list>
<data type="double"/>
<data type="double"/>
<data type="double"/>
</list>
</element>
<attribute name="reference_uvw">
<list>
<data type="double"/>
<data type="double"/>
<data type="double"/>
</list>
</attribute>
</choice>
</optional>
<optional>
<element name="mu">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="phi">
<ref name="distribution"/>
</element>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="energy">
<ref name="distribution"/>
</element>
</optional>
<optional>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<element name="write_initial">
<data type="boolean"/>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<attribute name="write_initial">
<data type="boolean"/>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<choice>
<element name="write_initial">
<data type="boolean"/>
</element>
<attribute name="write_initial">
<data type="boolean"/>
</attribute>
</choice>
</optional>
</interleave>
</optional>
</interleave>
</start>
<define name="distribution">
<interleave>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="interpolation">
<data type="string"/>
</element>
<attribute name="interpolation">
<data type="string"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</define>
</grammar>
</element>
</optional>
</zeroOrMore>
<optional>
<element name="state_point">
<choice>

2
src/relaxng/tallies.rnc
View file

@ -41,7 +41,7 @@ element tallies {
(element type { xsd:string } | attribute type { xsd:string }) &
(element threshold { xsd:double} | attribute threshold { xsd:double }) &
(element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )?
}?
}*
}* &
element assume_separate { xsd:boolean }?

8
src/relaxng/tallies.rng
View file

@ -145,10 +145,10 @@
<value>mesh</value>
<value>energy</value>
<value>energyout</value>
<value>delayedgroup</value>
<value>mu</value>
<value>polar</value>
<value>azimuthal</value>
<value>delayedgroup</value>
</choice>
</element>
<attribute name="type">
@ -162,10 +162,10 @@
<value>mesh</value>
<value>energy</value>
<value>energyout</value>
<value>delayedgroup</value>
<value>mu</value>
<value>polar</value>
<value>azimuthal</value>
<value>delayedgroup</value>
</choice>
</attribute>
</choice>
@ -208,7 +208,7 @@
</oneOrMore>
</list>
</element>
<optional>
<zeroOrMore>
<element name="trigger">
<interleave>
<choice>
@ -251,7 +251,7 @@
</optional>
</interleave>
</element>
</optional>
</zeroOrMore>
</interleave>
</element>
</zeroOrMore>

178
src/source.F90
View file

@ -2,6 +2,8 @@ module source
use bank_header, only: Bank
use constants
use distribution_univariate, only: Discrete
use distribution_multivariate, only: SpatialBox
use error, only: fatal_error
use geometry, only: find_cell
use geometry_header, only: BASE_UNIVERSE
@ -99,13 +101,9 @@ contains
type(Bank), intent(inout) :: site ! source site
integer :: i ! dummy loop index
integer :: n_source ! number of source distributions
real(8) :: c ! cumulative frequency
real(8) :: r(3) ! sampled coordinates
real(8) :: phi ! azimuthal angle
real(8) :: mu ! cosine of polar angle
real(8) :: p_min(3) ! minimum coordinates of source
real(8) :: p_max(3) ! maximum coordinates of source
real(8) :: a ! Arbitrary parameter 'a'
real(8) :: b ! Arbitrary parameter 'b'
logical :: found ! Does the source particle exist within geometry?
type(Particle) :: p ! Temporary particle for using find_cell
integer, save :: num_resamples = 0 ! Number of resamples encountered
@ -116,130 +114,80 @@ contains
! Set the random number generator to the source stream.
call prn_set_stream(STREAM_SOURCE)
! Sample position
select case (external_source%type_space)
case (SRC_SPACE_BOX)
! Sample from among multiple source distributions
n_source = size(external_source)
if (n_source > 1) then
r(1) = prn()*sum(external_source(:)%strength)
c = ZERO
do i = 1, n_source
c = c + external_source(i)%strength
if (r(1) < c) exit
end do
else
i = 1
end if
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Set particle defaults
call p%initialize()
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Coordinates sampled uniformly over a box
p_min = external_source%params_space(1:3)
p_max = external_source%params_space(4:6)
r = (/ (prn(), i = 1,3) /)
site%xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p%coord(1)%xyz = site%xyz
p%coord(1)%uvw = [ ONE, ZERO, ZERO ]
! Sample spatial distribution
site%xyz(:) = external_source(i)%space%sample()
! Now search to see if location exists in geometry
call find_cell(p, found)
if (.not. found) then
num_resamples = num_resamples + 1
if (num_resamples == MAX_EXTSRC_RESAMPLES) then
call fatal_error("Maximum number of external source spatial &
&resamples reached!")
! Fill p with needed data
p%coord(1)%xyz(:) = site%xyz
p%coord(1)%uvw(:) = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
if (.not. found) then
num_resamples = num_resamples + 1
if (num_resamples == MAX_EXTSRC_RESAMPLES) then
call fatal_error("Maximum number of external source spatial &
&resamples reached!")
end if
end if
! Check if spatial site is in fissionable material
select type (space => external_source(i)%space)
type is (SpatialBox)
if (space%only_fissionable) then
if (p%material == MATERIAL_VOID) then
found = .false.
elseif (.not. materials(p%material)%fissionable) then
found = .false.
end if
end if
end do
call p%clear()
end select
end do
case (SRC_SPACE_FISSION)
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not.found)
! Set particle defaults
call p%initialize()
! Coordinates sampled uniformly over a box
p_min = external_source%params_space(1:3)
p_max = external_source%params_space(4:6)
r = (/ (prn(), i = 1,3) /)
site%xyz = p_min + r*(p_max - p_min)
! Fill p with needed data
p%coord(1)%xyz = site%xyz
p%coord(1)%uvw = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
if (.not. found) then
num_resamples = num_resamples + 1
if (num_resamples == MAX_EXTSRC_RESAMPLES) then
call fatal_error("Maximum number of external source spatial &
&resamples reached!")
end if
cycle
end if
if (p%material == MATERIAL_VOID) then
found = .false.
cycle
end if
if (.not. materials(p%material)%fissionable) found = .false.
end do
call p%clear()
case (SRC_SPACE_POINT)
! Point source
site%xyz = external_source%params_space
end select
call p%clear()
! Sample angle
select case (external_source%type_angle)
case (SRC_ANGLE_ISOTROPIC)
! Sample isotropic distribution
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
site%uvw(1) = mu
site%uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site%uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
site%uvw(:) = external_source(i)%angle%sample()
case (SRC_ANGLE_MONO)
! Monodirectional source
site%uvw = external_source%params_angle
case default
call fatal_error("No angle distribution specified for external source!")
end select
! Sample energy distribution
select case (external_source%type_energy)
case (SRC_ENERGY_MONO)
! Monoenergtic source
site%E = external_source%params_energy(1)
if (site%E >= energy_max_neutron) then
! Check for monoenergetic source above maximum neutron energy
select type (energy => external_source(i)%energy)
type is (Discrete)
if (any(energy%x >= energy_max_neutron)) then
call fatal_error("Source energy above range of energies of at least &
&one cross section table")
end if
case (SRC_ENERGY_MAXWELL)
a = external_source%params_energy(1)
do
! Sample Maxwellian fission spectrum
site%E = maxwell_spectrum(a)
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
end do
case (SRC_ENERGY_WATT)
a = external_source%params_energy(1)
b = external_source%params_energy(2)
do
! Sample Watt fission spectrum
site%E = watt_spectrum(a, b)
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
end do
case default
call fatal_error("No energy distribution specified for external source!")
end select
do
! Sample energy spectrum
site%E = external_source(i)%energy%sample()
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
end do
! Set delayed group
site%delayed_group = 0
! Set the random number generator back to the tracking stream.
call prn_set_stream(STREAM_TRACKING)

21
src/source_header.F90
View file

@ -1,19 +1,20 @@
module source_header
use distribution_univariate, only: Distribution
use distribution_multivariate, only: UnitSphereDistribution, SpatialDistribution
implicit none
!===============================================================================
! EXTSOURCE describes an external source of neutrons for a fixed-source problem
! or for the starting source in a k eigenvalue problem
! SOURCEDISTRIBUTION describes an external source of particles for a
! fixed-source problem or for the starting source in a k eigenvalue problem
!===============================================================================
type ExtSource
integer :: type_space ! spacial distribution, e.g. 'box' or 'point'
integer :: type_angle ! angle distribution, e.g. 'isotropic'
integer :: type_energy ! energy distribution, e.g. 'Watt'
real(8), allocatable :: params_space(:) ! parameters for spatial distribution
real(8), allocatable :: params_angle(:) ! parameters for angle distribution
real(8), allocatable :: params_energy(:) ! parameters for energy distribution
end type ExtSource
type SourceDistribution
real(8) :: strength ! source strength
class(SpatialDistribution), allocatable :: space ! spatial distribution
class(UnitSphereDistribution), allocatable :: angle ! angle distribution
class(Distribution), allocatable :: energy ! energy distribution
end type SourceDistribution
end module source_header

2
src/state_point.F90
View file

@ -22,6 +22,7 @@ module state_point
use tally_header, only: TallyObject
use mesh_header, only: RegularMesh
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
use random_lcg, only: seed
#ifdef MPI
use message_passing
@ -268,7 +269,6 @@ contains
call write_dataset(filter_group, "type", "delayedgroup")
end select
call write_dataset(filter_group, "offset", tally%filters(j)%offset)
call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins)
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT .or. &

5
src/string.F90
View file

@ -3,7 +3,6 @@ module string
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL, &
OP_LEFT_PAREN, OP_RIGHT_PAREN, OP_COMPLEMENT, OP_INTERSECTION, OP_UNION
use error, only: fatal_error, warning
use global, only: master
use stl_vector, only: VectorInt
implicit none
@ -50,8 +49,8 @@ contains
if (i_end > 0) then
n = n + 1
if (i_end - i_start + 1 > len(words(n))) then
if (master) call warning("The word '" // string(i_start:i_end) &
&// "' is longer than the space allocated for it.")
call warning("The word '" // string(i_start:i_end) &
// "' is longer than the space allocated for it.")
end if
words(n) = string(i_start:i_end)
! reset indices

73
src/summary.F90
View file

@ -107,6 +107,7 @@ contains
integer :: i, j, k, m
integer, allocatable :: lattice_universes(:,:,:)
integer, allocatable :: cell_materials(:)
integer(HID_T) :: geom_group
integer(HID_T) :: cells_group, cell_group
integer(HID_T) :: surfaces_group, surface_group
@ -150,10 +151,24 @@ contains
select case (c%type)
case (CELL_NORMAL)
call write_dataset(cell_group, "fill_type", "normal")
if (c%material == MATERIAL_VOID) then
call write_dataset(cell_group, "material", -1)
if (size(c % material) == 1) then
if (c % material(1) == MATERIAL_VOID) then
call write_dataset(cell_group, "material", MATERIAL_VOID)
else
call write_dataset(cell_group, "material", &
materials(c % material(1)) % id)
end if
else
call write_dataset(cell_group, "material", materials(c%material)%id)
allocate(cell_materials(size(c % material)))
do j = 1, size(c % material)
if (c % material(j) == MATERIAL_VOID) then
cell_materials(j) = MATERIAL_VOID
else
cell_materials(j) = materials(c % material(j)) % id
end if
end do
call write_dataset(cell_group, "material", cell_materials)
deallocate(cell_materials)
end if
case (CELL_FILL)
@ -196,6 +211,8 @@ contains
end do
call write_dataset(cell_group, "region", adjustl(region_spec))
call write_dataset(cell_group, "distribcell_index", c % distribcell_index)
call close_group(cell_group)
end do CELL_LOOP
@ -355,16 +372,22 @@ contains
call write_dataset(lattice_group, "type", "rectangular")
! Write lattice dimensions, lower left corner, and pitch
call write_dataset(lattice_group, "dimension", lat%n_cells)
call write_dataset(lattice_group, "lower_left", lat%lower_left)
if (lat % is_3d) then
call write_dataset(lattice_group, "dimension", lat % n_cells)
call write_dataset(lattice_group, "lower_left", lat % lower_left)
else
call write_dataset(lattice_group, "dimension", lat % n_cells(1:2))
call write_dataset(lattice_group, "lower_left", lat % lower_left)
end if
! Write lattice universes.
allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), &
&lat%n_cells(3)))
do j = 1, lat%n_cells(1)
do k = 1, lat%n_cells(2)
do k = 0, lat%n_cells(2) - 1
do m = 1, lat%n_cells(3)
lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id
lattice_universes(j, k+1, m) = &
universes(lat%universes(j, lat%n_cells(2) - k, m))%id
end do
end do
end do
@ -484,8 +507,10 @@ contains
subroutine write_tallies(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i, j
integer :: i, j, k
integer :: i_list, i_xs
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer(HID_T) :: tallies_group
integer(HID_T) :: mesh_group
integer(HID_T) :: tally_group
@ -540,7 +565,6 @@ contains
filter_group = create_group(tally_group, "filter " // trim(to_str(j)))
! Write number of bins for this filter
call write_dataset(filter_group, "offset", t%filters(j)%offset)
call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins)
! Write filter bins
@ -664,6 +688,37 @@ contains
deallocate(str_array)
! Write explicit moment order strings for each score bin
k = 1
allocate(str_array(t%n_score_bins))
MOMENT_LOOP: do j = 1, t%n_user_score_bins
select case(t%score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
str_array(k) = 'P' // trim(to_str(t%moment_order(k)))
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, t%moment_order(k)
str_array(k) = 'P' // trim(to_str(n_order))
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
do n_order = 0, t%moment_order(k)
do nm_order = -n_order, n_order
str_array(k) = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
k = k + 1
end do
end do
case default
str_array(k) = ''
k = k + 1
end select
end do MOMENT_LOOP
call write_dataset(tally_group, "moment_orders", str_array)
deallocate(str_array)
call close_group(tally_group)
end do TALLY_METADATA

7
src/tally.F90
View file

@ -1871,6 +1871,7 @@ contains
integer :: j
integer :: n ! number of bins for single filter
integer :: offset ! offset for distribcell
integer :: distribcell_index ! index in distribcell arrays
real(8) :: E ! particle energy
real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively
type(TallyObject), pointer :: t
@ -1915,12 +1916,14 @@ contains
case (FILTER_DISTRIBCELL)
! determine next distribcell bin
distribcell_index = cells(t % filters(i) % int_bins(1)) &
% distribcell_index
matching_bins(i) = NO_BIN_FOUND
offset = 0
do j = 1, p % n_coord
if (cells(p % coord(j) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(j) % cell) % &
offset(t % filters(i) % offset)
offset(distribcell_index)
elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(j + 1) % lattice) % obj &
% are_valid_indices([&
@ -1928,7 +1931,7 @@ contains
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z])) then
offset = offset + lattices(p % coord(j + 1) % lattice) % obj % &
offset(t % filters(i) % offset, &
offset(distribcell_index, &
p % coord(j + 1) % lattice_x, &
p % coord(j + 1) % lattice_y, &
p % coord(j + 1) % lattice_z)

1
src/tally_header.F90
View file

@ -55,7 +55,6 @@ module tally_header
type TallyFilter
integer :: type = NONE
integer :: n_bins = 0
integer :: offset = 0 ! Only used for distribcell filters
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
end type TallyFilter

5
tests/input_set.py
View file

@ -1,4 +1,6 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
class InputSet(object):
@ -558,7 +560,8 @@ class InputSet(object):
self.settings.batches = 10
self.settings.inactive = 5
self.settings.particles = 100
self.settings.set_source_space('box', (-160, -160, -183, 160, 160, 183))
self.settings.source = Source(space=Box(
[-160, -160, -183], [160, 160, 183]))
def build_defualt_plots(self):
plot = openmc.Plot()

1
tests/test_distribmat/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
401b8be1b296db7f21ccae089c7ac480044d953b7264ca0ae8e34bb79e24cbb57195bcb568deda6f2f7e07366bbfac408a92306351b9169edd04499723707e1b

11
tests/test_distribmat/results_true.dat Normal file
View file

@ -0,0 +1,11 @@
k-combined:
1.309285E+00 1.263629E-02
Cell
ID = 11
Name =
Material = [2, 3, void, 2]
Region = -10000
Rotation = None
Translation = None
Offset = None
Distribcell index= 1

136
tests/test_distribmat/test_distribmat.py Normal file
View file

@ -0,0 +1,136 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
from openmc.stats import Box
from openmc.source import Source
class DistribmatTestHarness(PyAPITestHarness):
def _build_inputs(self):
####################
# Materials
####################
moderator = openmc.Material(material_id=1)
moderator.set_density('g/cc', 1.0)
moderator.add_nuclide('H-1', 2.0)
moderator.add_nuclide('O-16', 1.0)
dense_fuel = openmc.Material(material_id=2)
dense_fuel.set_density('g/cc', 4.5)
dense_fuel.add_nuclide('U-235', 1.0)
light_fuel = openmc.Material(material_id=3)
light_fuel.set_density('g/cc', 2.0)
light_fuel.add_nuclide('U-235', 1.0)
mats_file = openmc.MaterialsFile()
mats_file.default_xs = '71c'
mats_file.add_materials([moderator, dense_fuel, light_fuel])
mats_file.export_to_xml()
####################
# Geometry
####################
c1 = openmc.Cell(cell_id=1)
c1.fill = moderator
mod_univ = openmc.Universe(universe_id=1)
mod_univ.add_cell(c1)
r0 = openmc.ZCylinder(R=0.3)
c11 = openmc.Cell(cell_id=11)
c11.region = -r0
c11.fill = [dense_fuel, light_fuel, 'void', dense_fuel]
c12 = openmc.Cell(cell_id=12)
c12.region = +r0
c12.fill = moderator
fuel_univ = openmc.Universe(universe_id=11)
fuel_univ.add_cells((c11, c12))
lat = openmc.RectLattice(lattice_id=101)
lat.dimension = [2, 2]
lat.lower_left = [-2.0, -2.0]
lat.pitch = [2.0, 2.0]
lat.universes = [[fuel_univ]*2]*2
lat.outer = mod_univ
x0 = openmc.XPlane(x0=-3.0)
x1 = openmc.XPlane(x0=3.0)
y0 = openmc.YPlane(y0=-3.0)
y1 = openmc.YPlane(y0=3.0)
for s in [x0, x1, y0, y1]:
s.boundary_type = 'reflective'
c101 = openmc.Cell(cell_id=101)
c101.region = +x0 & -x1 & +y0 & -y1
c101.fill = lat
root_univ = openmc.Universe(universe_id=0)
root_univ.add_cell(c101)
geometry = openmc.Geometry()
geometry.root_universe = root_univ
geo_file = openmc.GeometryFile()
geo_file.geometry = geometry
geo_file.export_to_xml()
####################
# Settings
####################
sets_file = openmc.SettingsFile()
sets_file.batches = 5
sets_file.inactive = 0
sets_file.particles = 1000
sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1]))
sets_file.output = {'summary': True}
sets_file.export_to_xml()
####################
# Plots
####################
plots_file = openmc.PlotsFile()
plot = openmc.Plot(plot_id=1)
plot.basis = 'xy'
plot.color = 'cell'
plot.filename = 'cellplot'
plot.origin = (0, 0, 0)
plot.width = (7, 7)
plot.pixels = (400, 400)
plots_file.add_plot(plot)
plot = openmc.Plot(plot_id=2)
plot.basis = 'xy'
plot.color = 'mat'
plot.filename = 'matplot'
plot.origin = (0, 0, 0)
plot.width = (7, 7)
plot.pixels = (400, 400)
plots_file.add_plot(plot)
plots_file.export_to_xml()
def _get_results(self):
outstr = super(DistribmatTestHarness, self)._get_results()
su = openmc.Summary('summary.h5')
outstr += str(su.get_cell_by_id(11))
return outstr
def _cleanup(self):
f = os.path.join(os.getcwd(), 'plots.xml')
if os.path.exists(f):
os.remove(f)
super(DistribmatTestHarness, self)._cleanup()
if __name__ == '__main__':
harness = DistribmatTestHarness('statepoint.5.*')
harness.main()

2
tests/test_filter_azimuthal/inputs_true.dat
View file

@ -1 +1 @@
1d5f81d12f607f4a8436dfb65167e2a2be55dbf86fbc2cc465cec274671be5eaff517d781a4d40e264bb695e3c66c7ff61a650217c99de2ca8c15ca747fe6b80
57d6fd9cb5180c38efd2729a5dea0708cbd5fd0bf7dcf0c9d5c9cef5d818aeab5a926d03e70dedcf1b60d5740938fb3ba80e6ccdb09c661d159c0893da3bd593

2
tests/test_filter_cell/inputs_true.dat
View file

@ -1 +1 @@
caae173f01f7073d634a68a5c4ce97177423e13596a863976e1c40303dc8c05afed457d5a1aa0ae73627ec953143e4f9c1f45bdbd3b0cca76433062467d59777
f8359184c02fbab5dca5368689a84924066ab1fb09cae575588ceddd696d5461db577498df9959365d89fe933e9b338390e44e362c603c6f2aa5bcf4acc14b20

2
tests/test_filter_cellborn/inputs_true.dat
View file

@ -1 +1 @@
2f24eb86cda981982a8db5bb110c72e9cef542c06e15b748c1c7e459f96b0d8ba0978b51dffc006e813cd2e2ae1fa0357336f8322ae263189841afde01f0327b
8ae662f8881ce8cdec550069c6233c2c91e9a10f7200af6892cf6f2d77712ccfa17895dbd2eee02e6daf3d665c6ed84b29e17d89ff519e70c37b36d75a431d53

2
tests/test_filter_delayedgroup/inputs_true.dat
View file

@ -1 +1 @@
e771470681d3b4af57a70d148f5eb728df57f1fd7bdb5d6f76ac556b69f10bf0cdd5645fe488f1e17f07cbb72e9fb34af2fd7ede95c8e39c5ffa6ea9b6c5810e
a7c8ce7ffbc3a7b965d8a3077a4d9132130561afef19047b279b2d23198e248b09664856a092a32394894e19fef7708cebad99b3839d735c4e98ae0c9af58cb7

5
tests/test_filter_distribcell/test_filter_distribcell.py
View file

@ -10,10 +10,7 @@ from testing_harness import *
class DistribcellTestHarness(TestHarness):
def __init__(self):
self._sp_name = None
self._tallies = True
self._opts = None
self._args = None
super(DistribcellTestHarness, self).__init__(None, True)
def execute_test(self):
"""Run OpenMC with the appropriate arguments and check the outputs."""

2
tests/test_filter_energy/inputs_true.dat
View file

@ -1 +1 @@
49835200052ee1a4c9583a7bb4e9430de7d01a6d7a8d4f63ae37be9bbac4fd8c7ed9135ecbc4ab4cb075fd4d77b322265783463dd07c127decceee8c9fe25bfd
51d3e2c43f36712a7b26c5fa26e0e2ca6fb9af205af04f0f8cd44c6b100e36382417c2c63d711e4677ce3c1958d15072727d5fd32424a3f6eb08d1f3b1c7db5a

2
tests/test_filter_energyout/inputs_true.dat
View file

@ -1 +1 @@
183b4a06cbd0930cfa4f28d0c385cf3ae93ab97c171580c2ba9eec0e4b716102ad83f510d62258ef6769c446e72e1d5ba9f630b171ee239ec04c9bbd1e56742e
f0810606c5f947a9fe03bcfc87de3883ce46f59d8603e02ed30f853ebf301b2dc6bdcd109889801ada9e6e0b7be4932efeca97d4beea875af8c8e3ecb7511444

2
tests/test_filter_group_transfer/inputs_true.dat
View file

@ -1 +1 @@
461a6a4ec3b0b6dc7199c09fa8f527e51cc7a1ea4281a708f8b7bcdbb12f7162027928dfcef68a24eaf6c06037a68e151df9aac9eac6ce56617466f2f5270b71
c4d4334d44956d6dc9abe854a5e9403d7f8a87ffb04a15a3d128e8d18eb4111f46ca277b751e1b0e836d69527502f9abba115a4b2fc64c38da63a9d57968d860

2
tests/test_filter_material/inputs_true.dat
View file

@ -1 +1 @@
2fbd0986abff08126680925284929cf67bdad0cd564775197b78065ce3b0e6ad8094f5ca4d14ea69347db69a6cdcc9796a095178ae9b14a927b101f32f3cd0ec
7689b2c88391128377b7f9bfcda347a42f77d69d194186629fa965ecd3fc51be0bfd1ac92fb9d7551128d8b6ed5241ead4fb94b27ae29d80230863e78fbbcb68

2
tests/test_filter_mu/inputs_true.dat
View file

@ -1 +1 @@
cd2aeb24baafe9a904e697955990f6cffb5f25618fdf8c972775715bfe6e92bc259e36fd2b5addff8181439de58ad6f530972ec391e827f46146a2aaace34358
ecc649936e2cc364b079944f47e18fb81ec7290017b4bd5837e5aa1e24e1146df77897f44c7c2a88500e3f525566b51777cd9b84ec6a636f5883e411e4c1f75c

2
tests/test_filter_polar/inputs_true.dat
View file

@ -1 +1 @@
2de29e0a083af0722039ebc246469f26e260d604ab8b76314b2ac472bbd23004e031aec7afe3dbfc4c6112428846cd0cf0c1430e878832b9dab36463d026dee6
301824991a022884215609f39797a61933faf7ccacf81ad6bb883af08857563e8bd74ab946fc4fd072860168d77f76d0c76d1467375158072dce431fc6a1c449

2
tests/test_filter_universe/inputs_true.dat
View file

@ -1 +1 @@
51fcaf0aa527d1fd1022e5f312d4d25cd9bacc5fc9792d9f7702ee97cac43700a31f25be767a2cff769c37b5e1cdf3f922467977a9958fa34561e2a102bf8537
164804414f48a818c93e197f2901ce6ae375d88071a03e89c920dbc4462e7a2c8d2c85acf6560fcd6eb3d7c0c53d3b426ab1cc4b7721266fe8adec3e7231149e

12
tests/test_many_scores/results_true.dat
View file

@ -11,18 +11,6 @@ tally 1:
2.483728E+01
5.102293E-01
8.710841E-02
8.628000E+00
2.481430E+01
9.329009E-01
2.902534E-01
5.102293E-01
8.710841E-02
5.087118E-01
8.657086E-02
8.632000E+00
2.483728E+01
9.328366E-01
2.902108E-01
5.087118E-01
8.657086E-02
9.212024E+00

6
tests/test_many_scores/tallies.xml
View file

@ -4,9 +4,9 @@
<tally id="1">
<filter type="cell" bins="21" />
<scores>
flux total scatter nu-scatter scatter-2 scatter-p2 nu-scatter-2
nu-scatter-p2 transport n1n absorption nu-fission kappa-fission
flux-y2 total-y2 scatter-y2 nu-scatter-y2 events delayed-nu-fission
flux total scatter nu-scatter scatter-2 nu-scatter-2 transport n1n
absorption nu-fission kappa-fission flux-y2 total-y2 scatter-y2
nu-scatter-y2 events delayed-nu-fission
</scores>
</tally>

2
tests/test_mgxs_library_condense/inputs_true.dat
View file

@ -1 +1 @@
35f99f1973b3bf3efcec6c2dddf56d6679a15dab8582ab5336e86e4fdf90967ce91036e5c30c345decb994ab9133a906b82dc8fad0cdd3398a612d9aa05c1c77
53b1740921b71e4ead909ab9e4c25f7d43990fe7d7051fde6f66c39c0a6082177385640244010e1b9dbeaf5f34adf1627e9603088af729fadd6b589c19102edc

Some files were not shown because too many files have changed in this diff Show more