crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -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-04-18 14:48:52 -04:00
commit 88165bfdaa
184 changed files with 14327 additions and 14196 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -26,6 +26,7 @@ examples/python/**/*.xml
docs/build
docs/source/_images/*.pdf
docs/source/_images/*.aux
docs/source/pythonapi/generated/
# Source build
build

2
.travis.yml
View file

@ -27,7 +27,7 @@ before_install:
- conda config --set always_yes yes --set changeps1 no
- conda update -q conda
- conda info -a
- conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy h5py pandas
- conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy h5py=2.5 pandas
- source activate test-environment
# Install GCC, MPICH, HDF5, PHDF5

7
docs/source/_templates/myclass.rst Normal file
View file

@ -0,0 +1,7 @@
{{ fullname }}
{{ underline }}
.. currentmodule:: {{ module }}
.. autoclass:: {{ objname }}
:members:

40
docs/source/conf.py
View file

@ -24,13 +24,8 @@ 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)
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
# If extensions (or modules to document with autodoc) are in another directory,
@ -48,6 +43,8 @@ extensions = ['sphinx.ext.autodoc',
'sphinx.ext.napoleon',
'sphinx.ext.mathjax',
'sphinx.ext.autosummary',
'sphinx.ext.intersphinx',
'sphinx.ext.viewcode',
'sphinx_numfig',
'notebook_sphinxext']
@ -65,7 +62,7 @@ master_doc = 'index'
# General information about the project.
project = u'OpenMC'
copyright = u'2011-2015, Massachusetts Institute of Technology'
copyright = u'2011-2016, Massachusetts Institute of Technology'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
@ -122,20 +119,13 @@ pygments_style = 'tango'
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. Major themes that come with
# Sphinx are currently 'default' and 'sphinxdoc'.
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'
# The theme to use for HTML and HTML Help pages
if not on_rtd:
import sphinx_rtd_theme
html_theme = 'sphinx_rtd_theme'
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = ["_theme"]
html_logo = '_images/openmc200px.png'
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".
@ -248,4 +238,12 @@ latex_elements = {
#Autodocumentation Flags
#autodoc_member_order = "groupwise"
#autoclass_content = "both"
#autosummary_generate = []
autosummary_generate = True
napoleon_use_ivar = True
intersphinx_mapping = {
'python': ('https://docs.python.org/3', None),
'numpy': ('http://docs.scipy.org/doc/numpy/', None),
'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None)
}

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

@ -1,8 +0,0 @@
.. _pythonapi_ace:
==========
ACE Format
==========
.. automodule:: openmc.ace
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_cmfd:
====
CMFD
====
.. automodule:: openmc.cmfd
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_element:
=======
Element
=======
.. automodule:: openmc.element
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_energy_groups:
=============
Energy Groups
=============
.. automodule:: openmc.mgxs.groups
:members:

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

@ -141,15 +141,12 @@
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\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"
"import openmc.mgxs as mgxs"
]
},
{
@ -342,9 +339,11 @@
"settings_file.inactive = inactive\n",
"settings_file.particles = particles\n",
"settings_file.output = {'tallies': True}\n",
"\n",
"# Create an initial uniform spatial source distribution over fissionable zones\n",
"bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n",
"settings_file.source = Source(space=Box(\n",
" bounds[:3], bounds[3:], only_fissionable=True))\n",
"uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n",
"settings_file.source = openmc.source.Source(space=uniform_dist)\n",
"\n",
"# Export to \"settings.xml\"\n",
"settings_file.export_to_xml()"
@ -423,22 +422,24 @@
"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), ('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)])"
"\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",
")])"
]
},
"execution_count": 13,
@ -518,8 +519,8 @@
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.1\n",
" Git SHA1: 34381b40a9445a727e360873aaa6ef892af1cb6a\n",
" Date/Time: 2016-02-07 15:58:16\n",
" Git SHA1: eeb5091ca3a34cc85df73a3318cae2b6c7097413\n",
" Date/Time: 2016-04-13 11:24:09\n",
" MPI Processes: 1\n",
"\n",
" ===========================================================================\n",
@ -546,56 +547,56 @@
"\n",
" Bat./Gen. k Average k \n",
" ========= ======== ==================== \n",
" 1/1 1.19804 \n",
" 2/1 1.12945 \n",
" 3/1 1.15573 \n",
" 4/1 1.13929 \n",
" 5/1 1.16300 \n",
" 6/1 1.22117 \n",
" 7/1 1.19012 \n",
" 8/1 1.11299 \n",
" 9/1 1.16066 \n",
" 10/1 1.12566 \n",
" 11/1 1.20854 \n",
" 12/1 1.14691 1.17773 +/- 0.03082\n",
" 13/1 1.17204 1.17583 +/- 0.01789\n",
" 14/1 1.14148 1.16724 +/- 0.01529\n",
" 15/1 1.17272 1.16834 +/- 0.01189\n",
" 16/1 1.18575 1.17124 +/- 0.01014\n",
" 17/1 1.20498 1.17606 +/- 0.00983\n",
" 18/1 1.14754 1.17249 +/- 0.00923\n",
" 19/1 1.18141 1.17348 +/- 0.00820\n",
" 20/1 1.15074 1.17121 +/- 0.00768\n",
" 21/1 1.15914 1.17011 +/- 0.00703\n",
" 22/1 1.14586 1.16809 +/- 0.00673\n",
" 23/1 1.18999 1.16978 +/- 0.00642\n",
" 24/1 1.15101 1.16844 +/- 0.00609\n",
" 25/1 1.13791 1.16640 +/- 0.00602\n",
" 26/1 1.19791 1.16837 +/- 0.00597\n",
" 27/1 1.19818 1.17012 +/- 0.00587\n",
" 28/1 1.14160 1.16854 +/- 0.00576\n",
" 29/1 1.11487 1.16571 +/- 0.00614\n",
" 30/1 1.17538 1.16620 +/- 0.00584\n",
" 31/1 1.20210 1.16791 +/- 0.00581\n",
" 32/1 1.20078 1.16940 +/- 0.00574\n",
" 33/1 1.14624 1.16839 +/- 0.00558\n",
" 34/1 1.14618 1.16747 +/- 0.00542\n",
" 35/1 1.16866 1.16752 +/- 0.00520\n",
" 36/1 1.18565 1.16821 +/- 0.00504\n",
" 37/1 1.16824 1.16821 +/- 0.00485\n",
" 38/1 1.18299 1.16874 +/- 0.00471\n",
" 39/1 1.21418 1.17031 +/- 0.00480\n",
" 40/1 1.11167 1.16835 +/- 0.00504\n",
" 41/1 1.11545 1.16665 +/- 0.00516\n",
" 42/1 1.11114 1.16491 +/- 0.00529\n",
" 43/1 1.14227 1.16423 +/- 0.00517\n",
" 44/1 1.14104 1.16355 +/- 0.00506\n",
" 45/1 1.16756 1.16366 +/- 0.00492\n",
" 46/1 1.13065 1.16274 +/- 0.00487\n",
" 47/1 1.11251 1.16139 +/- 0.00492\n",
" 48/1 1.14731 1.16101 +/- 0.00481\n",
" 49/1 1.16691 1.16117 +/- 0.00469\n",
" 50/1 1.19679 1.16206 +/- 0.00465\n",
" 1/1 1.11184 \n",
" 2/1 1.15820 \n",
" 3/1 1.18468 \n",
" 4/1 1.17492 \n",
" 5/1 1.19645 \n",
" 6/1 1.18436 \n",
" 7/1 1.14070 \n",
" 8/1 1.15150 \n",
" 9/1 1.19202 \n",
" 10/1 1.17677 \n",
" 11/1 1.20272 \n",
" 12/1 1.21366 1.20819 +/- 0.00547\n",
" 13/1 1.15906 1.19181 +/- 0.01668\n",
" 14/1 1.14687 1.18058 +/- 0.01629\n",
" 15/1 1.14570 1.17360 +/- 0.01442\n",
" 16/1 1.13480 1.16713 +/- 0.01343\n",
" 17/1 1.17680 1.16852 +/- 0.01144\n",
" 18/1 1.16866 1.16853 +/- 0.00990\n",
" 19/1 1.19253 1.17120 +/- 0.00913\n",
" 20/1 1.18124 1.17220 +/- 0.00823\n",
" 21/1 1.19206 1.17401 +/- 0.00766\n",
" 22/1 1.17681 1.17424 +/- 0.00700\n",
" 23/1 1.17634 1.17440 +/- 0.00644\n",
" 24/1 1.13659 1.17170 +/- 0.00654\n",
" 25/1 1.17144 1.17169 +/- 0.00609\n",
" 26/1 1.20649 1.17386 +/- 0.00610\n",
" 27/1 1.11238 1.17024 +/- 0.00678\n",
" 28/1 1.18911 1.17129 +/- 0.00647\n",
" 29/1 1.14681 1.17000 +/- 0.00626\n",
" 30/1 1.12152 1.16758 +/- 0.00641\n",
" 31/1 1.12729 1.16566 +/- 0.00639\n",
" 32/1 1.15399 1.16513 +/- 0.00612\n",
" 33/1 1.13547 1.16384 +/- 0.00599\n",
" 34/1 1.17723 1.16440 +/- 0.00576\n",
" 35/1 1.09296 1.16154 +/- 0.00622\n",
" 36/1 1.19621 1.16287 +/- 0.00612\n",
" 37/1 1.12560 1.16149 +/- 0.00605\n",
" 38/1 1.17872 1.16211 +/- 0.00586\n",
" 39/1 1.17721 1.16263 +/- 0.00568\n",
" 40/1 1.13724 1.16178 +/- 0.00555\n",
" 41/1 1.18526 1.16254 +/- 0.00542\n",
" 42/1 1.13779 1.16177 +/- 0.00531\n",
" 43/1 1.15066 1.16143 +/- 0.00516\n",
" 44/1 1.12174 1.16026 +/- 0.00514\n",
" 45/1 1.17479 1.16068 +/- 0.00501\n",
" 46/1 1.14146 1.16014 +/- 0.00489\n",
" 47/1 1.20464 1.16135 +/- 0.00491\n",
" 48/1 1.15119 1.16108 +/- 0.00479\n",
" 49/1 1.17938 1.16155 +/- 0.00468\n",
" 50/1 1.15798 1.16146 +/- 0.00457\n",
" Creating state point statepoint.50.h5...\n",
"\n",
" ===========================================================================\n",
@ -605,27 +606,27 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 3.2100E-01 seconds\n",
" Reading cross sections = 7.4000E-02 seconds\n",
" Total time in simulation = 8.3830E+00 seconds\n",
" Time in transport only = 8.3670E+00 seconds\n",
" Time in inactive batches = 1.0330E+00 seconds\n",
" Time in active batches = 7.3500E+00 seconds\n",
" Time synchronizing fission bank = 4.0000E-03 seconds\n",
" Sampling source sites = 1.0000E-03 seconds\n",
" Total time for initialization = 4.6300E-01 seconds\n",
" Reading cross sections = 1.2100E-01 seconds\n",
" Total time in simulation = 1.6504E+01 seconds\n",
" Time in transport only = 1.6479E+01 seconds\n",
" Time in inactive batches = 1.9620E+00 seconds\n",
" Time in active batches = 1.4542E+01 seconds\n",
" Time synchronizing fission bank = 1.0000E-02 seconds\n",
" Sampling source sites = 4.0000E-03 seconds\n",
" SEND/RECV source sites = 3.0000E-03 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 8.7140E+00 seconds\n",
" Calculation Rate (inactive) = 24201.4 neutrons/second\n",
" Calculation Rate (active) = 13605.4 neutrons/second\n",
" Total time for finalization = 0.0000E+00 seconds\n",
" Total time elapsed = 1.6977E+01 seconds\n",
" Calculation Rate (inactive) = 12742.1 neutrons/second\n",
" Calculation Rate (active) = 6876.63 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.16131 +/- 0.00453\n",
" k-effective (Track-length) = 1.16206 +/- 0.00465\n",
" k-effective (Absorption) = 1.16096 +/- 0.00364\n",
" Combined k-effective = 1.16120 +/- 0.00325\n",
" k-effective (Collision) = 1.15984 +/- 0.00411\n",
" k-effective (Track-length) = 1.16146 +/- 0.00457\n",
" k-effective (Absorption) = 1.16177 +/- 0.00380\n",
" Combined k-effective = 1.16105 +/- 0.00364\n",
" Leakage Fraction = 0.00000 +/- 0.00000\n",
"\n"
]
@ -751,8 +752,8 @@
"\tDomain Type =\tcell\n",
"\tDomain ID =\t1\n",
"\tCross Sections [cm^-1]:\n",
" Group 1 [6.25e-07 - 20.0 MeV]:\t6.81e-01 +/- 1.88e-01%\n",
" Group 2 [0.0 - 6.25e-07 MeV]:\t1.40e+00 +/- 5.91e-01%\n",
" Group 1 [6.25e-07 - 20.0 MeV]:\t6.81e-01 +/- 2.69e-01%\n",
" Group 2 [0.0 - 6.25e-07 MeV]:\t1.40e+00 +/- 5.93e-01%\n",
"\n",
"\n",
"\n"
@ -780,7 +781,7 @@
{
"data": {
"text/html": [
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -795,19 +796,19 @@
" <tbody>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 1</td>\n",
" <td> 1</td>\n",
" <td> total</td>\n",
" <td> 0.668323</td>\n",
" <td> 0.001264</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>total</td>\n",
" <td>0.667787</td>\n",
" <td>0.001802</td>\n",
" </tr>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 1</td>\n",
" <td> 2</td>\n",
" <td> total</td>\n",
" <td> 1.293258</td>\n",
" <td> 0.007624</td>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>total</td>\n",
" <td>1.292013</td>\n",
" <td>0.007642</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -815,8 +816,8 @@
],
"text/plain": [
" cell group in nuclide mean std. dev.\n",
"1 1 1 total 0.668323 0.001264\n",
"0 1 2 total 1.293258 0.007624"
"1 1 1 total 0.667787 0.001802\n",
"0 1 2 total 1.292013 0.007642"
]
},
"execution_count": 20,
@ -891,7 +892,7 @@
{
"data": {
"text/html": [
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -908,23 +909,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 1</td>\n",
" <td> 0.000000</td>\n",
" <td> 0.000001</td>\n",
" <td> total</td>\n",
" <td> (((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td> 4.884981e-15</td>\n",
" <td> 0.011274</td>\n",
" <td>1</td>\n",
" <td>0.000000e+00</td>\n",
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>-3.774758e-15</td>\n",
" <td>0.011292</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 1</td>\n",
" <td> 0.000001</td>\n",
" <td> 20.000000</td>\n",
" <td> total</td>\n",
" <td> (((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td> 1.221245e-15</td>\n",
" <td> 0.001802</td>\n",
" <td>1</td>\n",
" <td>6.250000e-07</td>\n",
" <td>2.000000e+01</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>1.443290e-15</td>\n",
" <td>0.002570</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -935,9 +936,9 @@
"0 1 0.00e+00 6.25e-07 total \n",
"1 1 6.25e-07 2.00e+01 total \n",
"\n",
" score mean std. dev. \n",
"0 (((total / flux) - (absorption / flux)) - (sca... 4.88e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... 1.22e-15 1.80e-03 "
" score mean std. dev. \n",
"0 (((total / flux) - (absorption / flux)) - (sca... -3.77e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... 1.44e-15 2.57e-03 "
]
},
"execution_count": 23,
@ -970,7 +971,7 @@
{
"data": {
"text/html": [
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -987,23 +988,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 1</td>\n",
" <td> 0.000000</td>\n",
" <td> 0.000001</td>\n",
" <td> total</td>\n",
" <td> ((absorption / flux) / (total / flux))</td>\n",
" <td> 0.076219</td>\n",
" <td> 0.000651</td>\n",
" <td>1</td>\n",
" <td>0.000000e+00</td>\n",
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>((absorption / flux) / (total / flux))</td>\n",
" <td>0.076115</td>\n",
" <td>0.000649</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 1</td>\n",
" <td> 0.000001</td>\n",
" <td> 20.000000</td>\n",
" <td> total</td>\n",
" <td> ((absorption / flux) / (total / flux))</td>\n",
" <td> 0.019319</td>\n",
" <td> 0.000086</td>\n",
" <td>1</td>\n",
" <td>6.250000e-07</td>\n",
" <td>2.000000e+01</td>\n",
" <td>total</td>\n",
" <td>((absorption / flux) / (total / flux))</td>\n",
" <td>0.019263</td>\n",
" <td>0.000095</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1015,8 +1016,8 @@
"1 1 6.25e-07 2.00e+01 total \n",
"\n",
" score mean std. dev. \n",
"0 ((absorption / flux) / (total / flux)) 7.62e-02 6.51e-04 \n",
"1 ((absorption / flux) / (total / flux)) 1.93e-02 8.65e-05 "
"0 ((absorption / flux) / (total / flux)) 7.61e-02 6.49e-04 \n",
"1 ((absorption / flux) / (total / flux)) 1.93e-02 9.46e-05 "
]
},
"execution_count": 24,
@ -1042,7 +1043,7 @@
{
"data": {
"text/html": [
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1059,23 +1060,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 1</td>\n",
" <td> 0.000000</td>\n",
" <td> 0.000001</td>\n",
" <td> total</td>\n",
" <td> ((scatter / flux) / (total / flux))</td>\n",
" <td> 0.923781</td>\n",
" <td> 0.007714</td>\n",
" <td>1</td>\n",
" <td>0.000000e+00</td>\n",
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>((scatter / flux) / (total / flux))</td>\n",
" <td>0.923885</td>\n",
" <td>0.007736</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 1</td>\n",
" <td> 0.000001</td>\n",
" <td> 20.000000</td>\n",
" <td> total</td>\n",
" <td> ((scatter / flux) / (total / flux))</td>\n",
" <td> 0.980681</td>\n",
" <td> 0.002617</td>\n",
" <td>1</td>\n",
" <td>6.250000e-07</td>\n",
" <td>2.000000e+01</td>\n",
" <td>total</td>\n",
" <td>((scatter / flux) / (total / flux))</td>\n",
" <td>0.980737</td>\n",
" <td>0.003737</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1087,8 +1088,8 @@
"1 1 6.25e-07 2.00e+01 total \n",
"\n",
" score mean std. dev. \n",
"0 ((scatter / flux) / (total / flux)) 9.24e-01 7.71e-03 \n",
"1 ((scatter / flux) / (total / flux)) 9.81e-01 2.62e-03 "
"0 ((scatter / flux) / (total / flux)) 9.24e-01 7.74e-03 \n",
"1 ((scatter / flux) / (total / flux)) 9.81e-01 3.74e-03 "
]
},
"execution_count": 25,
@ -1121,7 +1122,7 @@
{
"data": {
"text/html": [
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<div>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1138,23 +1139,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td> 1</td>\n",
" <td> 0.000000</td>\n",
" <td> 0.000001</td>\n",
" <td> total</td>\n",
" <td> (((absorption / flux) / (total / flux)) + ((sc...</td>\n",
" <td> 1</td>\n",
" <td> 0.007741</td>\n",
" <td>1</td>\n",
" <td>0.000000e+00</td>\n",
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(((absorption / flux) / (total / flux)) + ((sc...</td>\n",
" <td>1</td>\n",
" <td>0.007763</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td> 1</td>\n",
" <td> 0.000001</td>\n",
" <td> 20.000000</td>\n",
" <td> total</td>\n",
" <td> (((absorption / flux) / (total / flux)) + ((sc...</td>\n",
" <td> 1</td>\n",
" <td> 0.002619</td>\n",
" <td>1</td>\n",
" <td>6.250000e-07</td>\n",
" <td>2.000000e+01</td>\n",
" <td>total</td>\n",
" <td>(((absorption / flux) / (total / flux)) + ((sc...</td>\n",
" <td>1</td>\n",
" <td>0.003739</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1166,8 +1167,8 @@
"1 1 6.25e-07 2.00e+01 total \n",
"\n",
" score mean std. dev. \n",
"0 (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00 7.74e-03 \n",
"1 (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00 2.62e-03 "
"0 (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00 7.76e-03 \n",
"1 (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00 3.74e-03 "
]
},
"execution_count": 26,
@ -1200,7 +1201,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
"version": "2.7.6"
}
},
"nbformat": 4,

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because it is too large Load diff

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

@ -1,8 +0,0 @@
.. _pythonapi_executor:
========
Executor
========
.. automodule:: openmc.executor
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_filter:
======
Filter
======
.. automodule:: openmc.filter
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_geometry:
========
Geometry
========
.. automodule:: openmc.geometry
:members:

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

@ -13,62 +13,267 @@ online. We recommend going through the modules from Codecademy_ and/or the
`Scipy lectures`_. The full API documentation serves to provide more information
on a given module or class.
**Handling nuclear data:**
------------------------------------
:mod:`openmc` -- Basic Functionality
------------------------------------
.. toctree::
:maxdepth: 1
Handling nuclear data
---------------------
ace
Classes
+++++++
**Creating input files:**
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
.. toctree::
:maxdepth: 1
openmc.XSdata
openmc.MGXSLibraryFile
cmfd
element
filter
geometry
material
mesh
nuclide
opencg_compatible
plots
settings
source
stats
surface
tallies
trigger
universe
Functions
+++++++++
**Running OpenMC:**
.. autosummary::
:toctree: generated
:nosignatures:
.. toctree::
:maxdepth: 1
openmc.ace.ascii_to_binary
executor
Simulation Settings
-------------------
**Post-processing:**
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
.. toctree::
:maxdepth: 1
openmc.Source
openmc.ResonanceScattering
openmc.SettingsFile
particle_restart
statepoint
summary
tallies
Material Specification
----------------------
**Multi-Group Cross Section Generation**
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
.. toctree::
:maxdepth: 1
openmc.Nuclide
openmc.Element
openmc.Macroscopic
openmc.Material
openmc.MaterialsFile
mgxs
energy_groups
mgxs_library
Building geometry
-----------------
**Example Jupyter Notebooks:**
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Plane
openmc.XPlane
openmc.YPlane
openmc.ZPlane
openmc.XCylinder
openmc.YCylinder
openmc.ZCylinder
openmc.Sphere
openmc.Cone
openmc.XCone
openmc.YCone
openmc.ZCone
openmc.Quadric
openmc.Halfspace
openmc.Intersection
openmc.Union
openmc.Complement
openmc.Cell
openmc.Universe
openmc.RectLattice
openmc.HexLattice
openmc.Geometry
openmc.GeometryFile
Many of the above classes are derived from several abstract classes:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Surface
openmc.Region
openmc.Lattice
Constructing Tallies
--------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Filter
openmc.Mesh
openmc.Trigger
openmc.Tally
openmc.TalliesFile
Coarse Mesh Finite Difference Acceleration
------------------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.CMFDMesh
openmc.CMFDFile
Plotting
--------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Plot
openmc.PlotsFile
Running OpenMC
--------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Executor
Post-processing
---------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.Particle
openmc.StatePoint
openmc.Summary
Various classes may be created when performing tally slicing and/or arithmetic:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.arithmetic.CrossScore
openmc.arithmetic.CrossNuclide
openmc.arithmetic.CrossFilter
openmc.arithmetic.AggregateScore
openmc.arithmetic.AggregateNuclide
openmc.arithmetic.AggregateFilter
---------------------------------
:mod:`openmc.stats` -- Statistics
---------------------------------
Univariate Probability Distributions
------------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.stats.Univariate
openmc.stats.Discrete
openmc.stats.Uniform
openmc.stats.Maxwell
openmc.stats.Watt
openmc.stats.Tabular
Angular Distributions
---------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.stats.UnitSphere
openmc.stats.PolarAzimuthal
openmc.stats.Isotropic
openmc.stats.Monodirectional
Spatial Distributions
---------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.stats.Spatial
openmc.stats.CartesianIndependent
openmc.stats.Box
openmc.stats.Point
----------------------------------------------------------
:mod:`openmc.mgxs` -- Multi-Group Cross Section Generation
----------------------------------------------------------
Energy Groups
-------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.mgxs.EnergyGroups
Multi-group Cross Sections
--------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.mgxs.MGXS
openmc.mgxs.AbsorptionXS
openmc.mgxs.CaptureXS
openmc.mgxs.Chi
openmc.mgxs.FissionXS
openmc.mgxs.NuFissionXS
openmc.mgxs.NuScatterXS
openmc.mgxs.NuScatterMatrixXS
openmc.mgxs.ScatterXS
openmc.mgxs.ScatterMatrixXS
openmc.mgxs.TotalXS
openmc.mgxs.TransportXS
Multi-group Cross Section Libraries
-----------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.mgxs.Library
-------------------------
Example Jupyter Notebooks
-------------------------
.. toctree::
:maxdepth: 1

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

@ -1,8 +0,0 @@
.. _pythonapi_material:
=========
Materials
=========
.. automodule:: openmc.material
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_mesh:
====
Mesh
====
.. automodule:: openmc.mesh
:members:

66
docs/source/pythonapi/mgxs.rst
View file

@ -1,66 +0,0 @@
.. _pythonapi_mgxs:
==========================
Multi-Group Cross Sections
==========================
.. currentmodule:: openmc.mgxs.mgxs
----------------------------
Summary of Available Classes
----------------------------
.. autosummary::
MGXS
AbsorptionXS
CaptureXS
Chi
FissionXS
NuFissionXS
NuScatterXS
NuScatterMatrixXS
ScatterXS
ScatterMatrixXS
TotalXS
TransportXS
-------------------
Class Documentation
-------------------
.. autoclass:: MGXS
:members:
.. autoclass:: AbsorptionXS
:members:
.. autoclass:: CaptureXS
:members:
.. autoclass:: Chi
:members:
.. autoclass:: FissionXS
:members:
.. autoclass:: NuFissionXS
:members:
.. autoclass:: NuScatterXS
:members:
.. autoclass:: NuScatterMatrixXS
:members:
.. autoclass:: ScatterXS
:members:
.. autoclass:: ScatterMatrixXS
:members:
.. autoclass:: TotalXS
:members:
.. autoclass:: TransportXS
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_mgxs_library:
============
MGXS Library
============
.. automodule:: openmc.mgxs.library
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_nuclide:
=======
Nuclide
=======
.. automodule:: openmc.nuclide
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_particle_restart:
================
Particle Restart
================
.. automodule:: openmc.particle_restart
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_plots:
=====
Plots
=====
.. automodule:: openmc.plots
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_settings:
========
Settings
========
.. automodule:: openmc.settings
:members:

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

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

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

@ -1,8 +0,0 @@
.. _pythonapi_statepoint:
==========
Statepoint
==========
.. automodule:: openmc.statepoint
:members:

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

@ -1,58 +0,0 @@
.. _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:

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

@ -1,8 +0,0 @@
.. _pythonapi_summary:
=======
Summary
=======
.. automodule:: openmc.summary
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_surface:
=======
Surface
=======
.. automodule:: openmc.surface
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_tallies:
=======
Tallies
=======
.. automodule:: openmc.tallies
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_trigger:
=======
Trigger
=======
.. automodule:: openmc.trigger
:members:

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

@ -1,8 +0,0 @@
.. _pythonapi_universe:
========
Universe
========
.. automodule:: openmc.universe
:members:

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

@ -1258,6 +1258,9 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
.. note:: The ``scattering`` attribute/sub-element is not used in the
multi-group :ref:`energy_mode`.
:element:
Specifies that a natural element is present in the material. The natural
@ -1293,6 +1296,9 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
.. note:: The ``scattering`` attribute/sub-element is not used in the
multi-group :ref:`energy_mode`.
:sab:
Associates an S(a,b) table with the material. This element has
attributes/sub-elements called ``name`` and ``xs``. The ``name`` attribute
@ -1301,6 +1307,8 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
:macroscopic:
The ``macroscopic`` element is similar to the ``nuclide`` element, but,
recognizes that some multi-group libraries may be providing material

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

@ -264,7 +264,7 @@ if run_mode == 'k-eigenvalue':
Accumulated sum and sum-of-squares for each global tally. The compound type
has fields named ``sum`` and ``sum_sq``.
**tallies_present** (*int*)
**/tallies_present** (*int*)
Flag indicated if tallies are present in the file.
@ -276,3 +276,69 @@ if (run_mode == 'k-eigenvalue' and source_present > 0)
``wgt``, ``xyz``, ``uvw``, ``E``, ``g``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.
**/runtime/total initialization** (*double*)
Time (in seconds on the master process) spent reading inputs, allocating
arrays, etc.
**/runtime/reading cross sections** (*double*)
Time (in seconds on the master process) spent loading cross section
libraries (this is a subset of initialization).
**/runtime/simulation** (*double*)
Time (in seconds on the master process) spent between initialization and
finalization.
**/runtime/transport** (*double*)
Time (in seconds on the master process) spent transporting particles.
**/runtime/inactive batches** (*double*)
Time (in seconds on the master process) spent in the inactive batches
(including non-transport activities like communcating sites).
**/runtime/active batches** (*double*)
Time (in seconds on the master process) spent in the active batches
(including non-transport activities like communicating sites).
**/runtime/synchronizing fission bank** (*double*)
Time (in seconds on the master process) spent sampling source particles
from fission sites and communicating them to other processes for load
balancing.
**/runtime/sampling source sites** (*double*)
Time (in seconds on the master process) spent sampling source particles
from fission sites.
**/runtime/SEND-RECV source sites** (*double*)
Time (in seconds on the master process) spent communicating source sites
between processes for load balancing.
**/runtime/accumulating tallies** (*double*)
Time (in seconds on the master process) spent communicating tally results
and evaluating their statistics.
**/runtime/CMFD** (*double*)
Time (in seconds on the master process) spent evaluating CMFD.
**/runtime/CMFD building matrices** (*double*)
Time (in seconds on the master process) spent buliding CMFD matrices.
**/runtime/CMFD solving matrices** (*double*)
Time (in seconds on the master process) spent solving CMFD matrices.
**/runtime/total** (*double*)
Total time spent (in seconds on the master process) in the program.

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

@ -293,6 +293,13 @@ The current revision of the summary file format is 1.
Filter offset (used for distribcell filter).
**/tallies/tally <uid>/filter <j>/paths** (*char[][]*)
The paths traversed through the CSG tree to reach each distribcell
instance (for 'distribcell' filters only). This consists of the integer
IDs for each universe, cell and lattice delimited by '->'. Each lattice
cell is specified by its (x,y) or (x,y,z) indices.
**/tallies/tally <uid>/filter <j>/n_bins** (*int*)
Number of bins for the j-th filter.

8
docs/source/usersguide/processing.rst
View file

@ -196,10 +196,10 @@ Data Extraction
A great deal of information is available in statepoint files (See
:ref:`usersguide_statepoint`), all of which is accessible through the Python
API. The ``openmc.statepoint`` module (see :ref:`pythonapi_statepoint`) provides
a class to load statepoints and access data as requested; it is used in many of
the provided plotting utilities, OpenMC's regression test suite, and can be used
in user-created scripts to carry out manipulations of the data.
API. The :class:`openmc.StatePoint` class can load statepoints and access data
as requested; it is used in many of the provided plotting utilities, OpenMC's
regression test suite, and can be used in user-created scripts to carry out
manipulations of the data.
An :ref:`example IPython notebook <notebook_post_processing>` demonstrates how
to extract data from a statepoint using the Python API.

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

@ -1,6 +1,5 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -94,7 +93,12 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box([-4, -4, -4], [4, 4, 4]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-4., -4., -4., 4., 4., 4.]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
@ -109,29 +113,20 @@ energyout_filter = openmc.Filter(type='energyout', bins=[0., 20.])
# Instantiate the first Tally
first_tally = openmc.Tally(tally_id=1, name='first tally')
first_tally.add_filter(cell_filter)
scores = ['total', 'scatter', 'nu-scatter', \
first_tally.filters = [cell_filter]
scores = ['total', 'scatter', 'nu-scatter',
'absorption', 'fission', 'nu-fission']
for score in scores:
first_tally.add_score(score)
first_tally.scores = scores
# Instantiate the second Tally
second_tally = openmc.Tally(tally_id=2, name='second tally')
second_tally.add_filter(cell_filter)
second_tally.add_filter(energy_filter)
scores = ['total', 'scatter', 'nu-scatter', \
'absorption', 'fission', 'nu-fission']
for score in scores:
second_tally.add_score(score)
second_tally.filters = [cell_filter, energy_filter]
second_tally.scores = scores
# Instantiate the third Tally
third_tally = openmc.Tally(tally_id=3, name='third tally')
third_tally.add_filter(cell_filter)
third_tally.add_filter(energy_filter)
third_tally.add_filter(energyout_filter)
scores = ['scatter', 'nu-scatter', 'nu-fission']
for score in scores:
third_tally.add_score(score)
third_tally.filters = [cell_filter, energy_filter, energyout_filter]
third_tally.scores = ['scatter', 'nu-scatter', 'nu-fission']
# Instantiate a TalliesFile, register all Tallies, and export to XML
tallies_file = openmc.TalliesFile()

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

@ -1,8 +1,5 @@
import numpy as np
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -119,7 +116,11 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(*outer_cube.bounding_box))
# Create an initial uniform spatial source distribution over fissionable zones
uniform_dist = openmc.stats.Box(*outer_cube.bounding_box, only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
###############################################################################

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

@ -1,6 +1,4 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -126,8 +124,12 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.keff_trigger = {'type' : 'std_dev', 'threshold' : 5E-4}
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100
@ -166,8 +168,8 @@ plot_file.export_to_xml()
# Instantiate a distribcell Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(openmc.Filter(type='distribcell', bins=[cell2.id]))
tally.add_score('total')
tally.filters = [openmc.Filter(type='distribcell', bins=[cell2.id])]
tally.scores = ['total']
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()

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

@ -1,6 +1,4 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -137,8 +135,12 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
@ -175,8 +177,8 @@ mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(mesh_filter)
tally.add_score('total')
tally.filters = [mesh_filter]
tally.scores = ['total']
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()

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

@ -1,6 +1,4 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -127,8 +125,12 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(
[-1, -1, -1], [1, 1, 1]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100
settings_file.export_to_xml()
@ -167,13 +169,13 @@ mesh_filter.mesh = mesh
# Instantiate tally Trigger
trigger = openmc.Trigger(trigger_type='rel_err', threshold=1E-2)
trigger.add_score('all')
trigger.scores = ['all']
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(mesh_filter)
tally.add_score('total')
tally.add_trigger(trigger)
tally.filters = [mesh_filter]
tally.scores = ['total']
tally.triggers = [trigger]
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()

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

@ -1,6 +1,4 @@
import openmc
from openmc.source import Source
from openmc.stats import Box
###############################################################################
# Simulation Input File Parameters
@ -170,8 +168,12 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(
[-0.62992, -0.62992, -1], [0.62992, 0.62992, 1]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
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]
@ -196,11 +198,8 @@ mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='tally 1')
tally.add_filter(energy_filter)
tally.add_filter(mesh_filter)
tally.add_score('flux')
tally.add_score('fission')
tally.add_score('nu-fission')
tally.filters = [energy_filter, mesh_filter]
tally.scores = ['flux', 'fission', 'nu-fission']
# Instantiate a TalliesFile, register all Tallies, and export to XML
tallies_file = openmc.TalliesFile()

12
examples/python/pincell_multigroup/build-xml.py
View file

@ -1,8 +1,6 @@
import numpy as np
import openmc
import openmc.mgxs
from openmc.source import Source
from openmc.stats import Box
import numpy as np
###############################################################################
# Simulation Input File Parameters
@ -145,7 +143,13 @@ settings_file.cross_sections = "./mg_cross_sections.xml"
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box([-0.63, -0.63, -1.], [0.63, 0.63, 1.]))
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.63, -0.63, -1, 0.63, 0.63, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()
###############################################################################
# Exporting to OpenMC tallies.xml File

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

@ -1,8 +1,5 @@
import numpy as np
import openmc
from openmc.stats import Box
from openmc.source import Source
###############################################################################
# Simulation Input File Parameters
@ -86,5 +83,10 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box(*cell.region.bounding_box))
# Create an initial uniform spatial source distribution over fissionable zones
uniform_dist = openmc.stats.Box(*cell.region.bounding_box,
only_fissionable=True)
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.export_to_xml()

5
openmc/__init__.py
View file

@ -1,3 +1,5 @@
from openmc.cell import *
from openmc.lattice import *
from openmc.element import *
from openmc.geometry import *
from openmc.nuclide import *
@ -17,6 +19,9 @@ from openmc.cmfd import *
from openmc.executor import *
from openmc.statepoint import *
from openmc.summary import *
from openmc.region import *
from openmc.source import *
from openmc.particle_restart import *
try:
from openmc.opencg_compatible import *

447
openmc/cell.py Normal file
View file

@ -0,0 +1,447 @@
from collections import OrderedDict, Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import warnings
import openmc
import openmc.checkvalue as cv
from openmc.surface import Halfspace
from openmc.region import Region, Intersection, Complement
if sys.version_info[0] >= 3:
basestring = str
# A static variable for auto-generated Cell IDs
AUTO_CELL_ID = 10000
def reset_auto_cell_id():
global AUTO_CELL_ID
AUTO_CELL_ID = 10000
class Cell(object):
"""A region of space defined as the intersection of half-space created by
quadric surfaces.
Parameters
----------
cell_id : int, optional
Unique identifier for the cell. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the cell. If not specified, the name is the empty string.
Attributes
----------
id : int
Unique identifier for the cell
name : str
Name of the cell
fill : openmc.Material or openmc.Universe or openmc.Lattice or 'void' or iterable of openmc.Material
Indicates what the region of space is filled with
region : openmc.Region
Region of space that is assigned to the cell.
rotation : numpy.ndarray
If the cell is filled with a universe, this array specifies the angles
in degrees about the x, y, and z axes that the filled universe should be
rotated.
translation : numpy.ndarray
If the cell is filled with a universe, this array specifies a vector
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
"""
def __init__(self, cell_id=None, name=''):
# Initialize Cell class attributes
self.id = cell_id
self.name = name
self._fill = None
self._type = None
self._region = None
self._rotation = None
self._translation = None
self._offsets = None
self._distribcell_index = None
def __eq__(self, other):
if not isinstance(other, Cell):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif self.fill != other.fill:
return False
elif self.region != other.region:
return False
elif self.rotation != other.rotation:
return False
elif self.translation != other.translation:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Cell\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
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, (openmc.Universe, openmc.Lattice)):
string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t',
self._fill._id)
else:
string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill)
string += '{0: <16}{1}{2}\n'.format('\tRegion', '=\t', self._region)
string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t',
self._rotation)
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
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def fill(self):
return self._fill
@property
def fill_type(self):
if isinstance(self.fill, openmc.Material):
return 'material'
elif isinstance(self.fill, openmc.Universe):
return 'universe'
elif isinstance(self.fill, openmc.Lattice):
return 'lattice'
else:
return None
@property
def region(self):
return self._region
@property
def rotation(self):
return self._rotation
@property
def translation(self):
return self._translation
@property
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:
global AUTO_CELL_ID
self._id = AUTO_CELL_ID
AUTO_CELL_ID += 1
else:
cv.check_type('cell ID', cell_id, Integral)
cv.check_greater_than('cell ID', cell_id, 0, equality=True)
self._id = cell_id
@name.setter
def name(self, name):
if name is not None:
cv.check_type('cell name', name, basestring)
self._name = name
else:
self._name = ''
@fill.setter
def fill(self, fill):
if isinstance(fill, basestring):
if fill.strip().lower() == 'void':
self._type = 'void'
else:
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
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, openmc.Universe):
self._type = 'fill'
elif isinstance(fill, openmc.Lattice):
self._type = 'lattice'
else:
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
self._fill = fill
@rotation.setter
def rotation(self, rotation):
cv.check_type('cell rotation', rotation, Iterable, Real)
cv.check_length('cell rotation', rotation, 3)
self._rotation = rotation
@translation.setter
def translation(self, translation):
cv.check_type('cell translation', translation, Iterable, Real)
cv.check_length('cell translation', translation, 3)
self._translation = translation
@offsets.setter
def offsets(self, offsets):
cv.check_type('cell offsets', offsets, Iterable)
self._offsets = offsets
@region.setter
def region(self, region):
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.
.. deprecated:: 0.7.1
Use the :attr:`Cell.region` property to directly specify a Region
expression.
Parameters
----------
surface : openmc.Surface
Quadric surface dividing space
halfspace : {-1, 1}
Indicate whether the negative or positive half-space is to be used
"""
warnings.warn("Cell.add_surface(...) has been deprecated and may be "
"removed in a future version. The region for a Cell "
"should be defined using the region property directly.",
DeprecationWarning)
if not isinstance(surface, openmc.Surface):
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
if halfspace not in [-1, +1]:
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" with halfspace ' \
'"{2}" since it is not +/-1'.format(surface, self._id, halfspace)
raise ValueError(msg)
# If no region has been assigned, simply use the half-space. Otherwise,
# take the intersection of the current region and the half-space
# specified
region = +surface if halfspace == 1 else -surface
if self.region is None:
self.region = region
else:
if isinstance(self.region, Intersection):
self.region.nodes.append(region)
else:
self.region = Intersection(self.region, region)
def get_cell_instance(self, path, distribcell_index):
# If the Cell is filled by a Material
if self._type == 'normal' or self._type == 'void':
offset = 0
# If the Cell is filled by a Universe
elif self._type == 'fill':
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_cell_instance(path, distribcell_index)
return offset
def get_all_nuclides(self):
"""Return all nuclides contained in the cell
Returns
-------
nuclides : dict
Dictionary whose keys are nuclide names and values are 2-tuples of
(nuclide, density)
"""
nuclides = OrderedDict()
if self._type != 'void':
nuclides.update(self._fill.get_all_nuclides())
return nuclides
def get_all_cells(self):
"""Return all cells that are contained within this one if it is filled with a
universe or lattice
Returns
-------
cells : dict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
"""
cells = OrderedDict()
if self._type == 'fill' or self._type == 'lattice':
cells.update(self._fill.get_all_cells())
return cells
def get_all_materials(self):
"""Return all materials that are contained within the cell
Returns
-------
materials : dict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
"""
materials = OrderedDict()
if self.fill_type == 'material':
materials[self.fill.id] = self.fill
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells()
for cell_id, cell in cells.items():
materials.update(cell.get_all_materials())
return materials
def get_all_universes(self):
"""Return all universes that are contained within this one if any of
its cells are filled with a universe or lattice.
Returns
-------
universes : dict
Dictionary whose keys are universe IDs and values are
:class:`Universe` instances
"""
universes = OrderedDict()
if self._type == 'fill':
universes[self._fill._id] = self._fill
universes.update(self._fill.get_all_universes())
elif self._type == 'lattice':
universes.update(self._fill.get_all_universes())
return universes
def create_xml_subelement(self, xml_element):
element = ET.Element("cell")
element.set("id", str(self.id))
if len(self._name) > 0:
element.set("name", str(self.name))
if isinstance(self.fill, basestring):
element.set("material", "void")
elif 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, (openmc.Universe, openmc.Lattice)):
element.set("fill", str(self.fill.id))
self.fill.create_xml_subelement(xml_element)
else:
element.set("fill", str(self.fill))
self.fill.create_xml_subelement(xml_element)
if self.region is not None:
# Set the region attribute with the region specification
element.set("region", str(self.region))
# Only surfaces that appear in a region are added to the geometry
# file, so the appropriate check is performed here. First we create
# a function which is called recursively to navigate through the CSG
# tree. When it reaches a leaf (a Halfspace), it creates a <surface>
# element for the corresponding surface if none has been created
# thus far.
def create_surface_elements(node, element):
if isinstance(node, Halfspace):
path = './surface[@id=\'{0}\']'.format(node.surface.id)
if xml_element.find(path) is None:
surface_subelement = node.surface.create_xml_subelement()
xml_element.append(surface_subelement)
elif isinstance(node, Complement):
create_surface_elements(node.node, element)
else:
for subnode in node.nodes:
create_surface_elements(subnode, element)
# Call the recursive function from the top node
create_surface_elements(self.region, xml_element)
if self.translation is not None:
element.set("translation", ' '.join(map(str, self.translation)))
if self.rotation is not None:
element.set("rotation", ' '.join(map(str, self.rotation)))
return element

72
openmc/checkvalue.py
View file

@ -41,25 +41,36 @@ def check_type(name, value, expected_type, expected_iter_type=None):
Description of value being checked
value : object
Object to check type of
expected_type : type
expected_type : type or Iterable of type
type to check object against
expected_iter_type : type or None, optional
expected_iter_type : type or Iterable of type or None, optional
Expected type of each element in value, assuming it is iterable. If
None, no check will be performed.
"""
if not _isinstance(value, expected_type):
msg = 'Unable to set "{0}" to "{1}" which is not of type "{2}"'.format(
name, value, expected_type.__name__)
if isinstance(expected_type, Iterable):
msg = 'Unable to set "{0}" to "{1}" which is not one of the ' \
'following types: "{2}"'.format(name, value, ', '.join(
[t.__name__ for t in expected_type]))
else:
msg = 'Unable to set "{0}" to "{1}" which is not of type "{2}"'.format(
name, value, expected_type.__name__)
raise ValueError(msg)
if expected_iter_type:
for item in value:
if not _isinstance(item, expected_iter_type):
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
if isinstance(expected_iter_type, Iterable):
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'one of the following types: "{2}"'.format(
name, value, ', '.join([t.__name__ for t in
expected_iter_type]))
else:
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
raise ValueError(msg)
@ -245,3 +256,50 @@ def check_greater_than(name, value, minimum, equality=False):
msg = 'Unable to set "{0}" to "{1}" since it is less than ' \
'or equal to "{2}"'.format(name, value, minimum)
raise ValueError(msg)
class CheckedList(list):
"""A list for which each element is type-checked as it's added
Parameters
----------
expected_type : type or Iterable of type
Type(s) which each element should be
name : str
Name of data being checked
items : Iterable, optional
Items to initialize the list with
"""
def __init__(self, expected_type, name, items=[]):
self.expected_type = expected_type
self.name = name
for item in items:
self.append(item)
def append(self, item):
"""Append item to list
Parameters
----------
item : object
Item to append
"""
check_type(self.name, item, self.expected_type)
super(CheckedList, self).append(item)
def insert(self, index, item):
"""Insert item before index
Parameters
----------
index : int
Index in list
item : object
Item to insert
"""
check_type(self.name, item, self.expected_type)
super(CheckedList, self).insert(index, item)

4
openmc/cmfd.py
View file

@ -69,7 +69,7 @@ class CMFDMesh(object):
to any tallies far away from fission source neutron regions. A ``2``
must be used to identify any fission source region.
"""
"""
def __init__(self):
self._lower_left = None
@ -219,7 +219,7 @@ class CMFDFile(object):
inner tolerance for Gauss-Seidel iterations when performing CMFD.
ktol : float
Tolerance on the eigenvalue when performing CMFD power iteration
cmfd_mesh : CMFDMesh
cmfd_mesh : openmc.CMFDMesh
Structured mesh to be used for acceleration
norm : float
Normalization factor applied to the CMFD fission source distribution

2
openmc/element.py
View file

@ -24,7 +24,7 @@ class Element(object):
Chemical symbol of the element, e.g. Pu
xs : str
Cross section identifier, e.g. 71c
scattering : 'data' or 'iso-in-lab' or None
scattering : {'data', 'iso-in-lab', None}
The type of angular scattering distribution to use
"""

3
openmc/executor.py
View file

@ -27,7 +27,8 @@ class Executor(object):
# Launch a subprocess to run OpenMC
p = subprocess.Popen(command, shell=True,
cwd=self._working_directory,
stdout=subprocess.PIPE)
stdout=subprocess.PIPE,
universal_newlines=True)
# Capture and re-print OpenMC output in real-time
while True:

44
openmc/filter.py
View file

@ -40,11 +40,14 @@ class Filter(object):
The bins for the filter
num_bins : Integral
The number of filter bins
mesh : Mesh or None
mesh : openmc.Mesh or None
A Mesh object for 'mesh' type filters.
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
distribcell_paths : list of str
The paths traversed through the CSG tree to reach each distribcell
instance (for 'distribcell' filters only)
"""
@ -56,6 +59,7 @@ class Filter(object):
self._bins = None
self._mesh = None
self._stride = None
self._distribcell_paths = None
if type is not None:
self.type = type
@ -110,6 +114,7 @@ class Filter(object):
clone._num_bins = self.num_bins
clone._mesh = copy.deepcopy(self.mesh, memo)
clone._stride = self.stride
clone._distribcell_paths = copy.deepcopy(self.distribcell_paths)
memo[id(self)] = clone
@ -152,6 +157,10 @@ class Filter(object):
def stride(self):
return self._stride
@property
def distribcell_paths(self):
return self._distribcell_paths
@type.setter
def type(self, type):
if type is None:
@ -246,12 +255,17 @@ class Filter(object):
self._stride = stride
@distribcell_paths.setter
def distribcell_paths(self, distribcell_paths):
cv.check_iterable_type('distribcell_paths', distribcell_paths, str)
self._distribcell_paths = distribcell_paths
def can_merge(self, other):
"""Determine if filter can be merged with another.
Parameters
----------
other : Filter
other : openmc.Filter
Filter to compare with
Returns
@ -296,12 +310,12 @@ class Filter(object):
Parameters
----------
other : Filter
other : openmc.Filter
Filter to merge with
Returns
-------
merged_filter : Filter
merged_filter : openmc.Filter
Filter resulting from the merge
"""
@ -341,7 +355,7 @@ class Filter(object):
Parameters
----------
other : Filter
other : openmc.Filter
The filter to query as a subset of this filter
Returns
@ -505,8 +519,8 @@ class Filter(object):
"""Builds a Pandas DataFrame for the Filter's bins.
This method constructs a Pandas DataFrame object for the filter with
columns annotated by filter bin information. This is a helper method
for the Tally.get_pandas_dataframe(...) method.
columns annotated by filter bin information. This is a helper method for
:meth:`Tally.get_pandas_dataframe`.
This capability has been tested for Pandas >=0.13.1. However, it is
recommended to use v0.16 or newer versions of Pandas since this method
@ -516,7 +530,7 @@ class Filter(object):
----------
data_size : Integral
The total number of bins in the tally corresponding to this filter
summary : None or Summary
summary : None or openmc.Summary
An optional Summary object to be used to construct columns for
distribcell tally filters (default is None). The geometric
information in the Summary object is embedded into a Multi-index
@ -632,18 +646,10 @@ class Filter(object):
# offsets to OpenCG LocalCoords linked lists
offsets_to_coords = {}
# Use OpenCG to compute LocalCoords linked list for
# each region and store in dictionary
for region in range(num_regions):
for offset, path in enumerate(self.distribcell_paths):
region = opencg_geometry.get_region_from_path(path)
coords = opencg_geometry.find_region(region)
path = opencg.get_path(coords)
cell_id = path[-1]
# 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_cell_instance(path)
offsets_to_coords[offset] = coords
offsets_to_coords[offset] = coords
# Each distribcell offset is a DataFrame bin
# Unravel the paths into DataFrame columns

36
openmc/geometry.py
View file

@ -17,7 +17,7 @@ class Geometry(object):
Attributes
----------
root_universe : openmc.universe.Universe
root_universe : openmc.Universe
Root universe which contains all others
"""
@ -63,15 +63,19 @@ class Geometry(object):
"""
# Extract the cell id from the path
last_index = path.rfind('>')
cell_id = int(path[last_index+1:])
# Find the distribcell index of the cell.
cells = self.get_all_cells()
for cell in cells:
if cell.id == path[-1]:
if cell.id == cell_id:
distribcell_index = cell.distribcell_index
break
else:
raise RuntimeError('Could not find cell {} specified in a \
distribcell filter'.format(path[-1]))
distribcell filter'.format(cell_id))
# Return memoize'd offset if possible
if (path, distribcell_index) in self._offsets:
@ -91,7 +95,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Cell
list of openmc.Cell
Cells in the geometry
"""
@ -112,7 +116,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Universe
list of openmc.Universe
Universes in the geometry
"""
@ -132,7 +136,7 @@ class Geometry(object):
Returns
-------
list of openmc.nuclide.Nuclide
list of openmc.Nuclide
Nuclides in the geometry
"""
@ -150,7 +154,7 @@ class Geometry(object):
Returns
-------
list of openmc.material.Material
list of openmc.Material
Materials in the geometry
"""
@ -173,7 +177,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Cell
list of openmc.Cell
Cells filled by Materials in the geometry
"""
@ -194,7 +198,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Universe
list of openmc.Universe
Universes with non-fill cells
"""
@ -217,7 +221,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Lattice
list of openmc.Lattice
Lattices in the geometry
"""
@ -248,7 +252,7 @@ class Geometry(object):
Returns
-------
list of openmc.material.Material
list of openmc.Material
Materials matching the queried name
"""
@ -288,7 +292,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Cell
list of openmc.Cell
Cells matching the queried name
"""
@ -328,7 +332,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Cell
list of openmc.Cell
Cells with fills matching the queried name
"""
@ -368,7 +372,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Universe
list of openmc.Universe
Universes matching the queried name
"""
@ -408,7 +412,7 @@ class Geometry(object):
Returns
-------
list of openmc.universe.Lattice
list of openmc.Lattice
Lattices matching the queried name
"""
@ -440,7 +444,7 @@ class GeometryFile(object):
Attributes
----------
geometry : Geometry
geometry : openmc.Geometry
The geometry to be used
"""

871
openmc/lattice.py Normal file
View file

@ -0,0 +1,871 @@
import abc
from collections import OrderedDict, Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import numpy as np
import openmc.checkvalue as cv
from openmc.universe import Universe, AUTO_UNIVERSE_ID
if sys.version_info[0] >= 3:
basestring = str
class Lattice(object):
"""A repeating structure wherein each element is a universe.
Parameters
----------
lattice_id : int, optional
Unique identifier for the lattice. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the lattice. If not specified, the name is the empty string.
Attributes
----------
id : int
Unique identifier for the lattice
name : str
Name of the lattice
pitch : float
Pitch of the lattice in cm
outer : int
The unique identifier of a universe to fill all space outside the
lattice
universes : numpy.ndarray of openmc.Universe
An array of universes filling each element of the lattice
"""
# This is an abstract class which cannot be instantiated
__metaclass__ = abc.ABCMeta
def __init__(self, lattice_id=None, name=''):
# Initialize Lattice class attributes
self.id = lattice_id
self.name = name
self._pitch = None
self._outer = None
self._universes = None
def __eq__(self, other):
if not isinstance(other, Lattice):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif self.pitch != other.pitch:
return False
elif self.outer != other.outer:
return False
elif self.universes != other.universes:
return False
else:
return True
def __ne__(self, other):
return not self == other
@property
def id(self):
return self._id
@property
def name(self):
return self._name
@property
def pitch(self):
return self._pitch
@property
def outer(self):
return self._outer
@property
def universes(self):
return self._universes
@id.setter
def id(self, lattice_id):
if lattice_id is None:
global AUTO_UNIVERSE_ID
self._id = AUTO_UNIVERSE_ID
AUTO_UNIVERSE_ID += 1
else:
cv.check_type('lattice ID', lattice_id, Integral)
cv.check_greater_than('lattice ID', lattice_id, 0, equality=True)
self._id = lattice_id
@name.setter
def name(self, name):
if name is not None:
cv.check_type('lattice name', name, basestring)
self._name = name
else:
self._name = ''
@outer.setter
def outer(self, outer):
cv.check_type('outer universe', outer, Universe)
self._outer = outer
@universes.setter
def universes(self, universes):
cv.check_iterable_type('lattice universes', universes, Universe,
min_depth=2, max_depth=3)
self._universes = np.asarray(universes)
def get_unique_universes(self):
"""Determine all unique universes in the lattice
Returns
-------
universes : collections.OrderedDict
Dictionary whose keys are universe IDs and values are
:class:`Universe` instances
"""
univs = OrderedDict()
for k in range(len(self._universes)):
for j in range(len(self._universes[k])):
if isinstance(self._universes[k][j], Universe):
u = self._universes[k][j]
univs[u._id] = u
else:
for i in range(len(self._universes[k][j])):
u = self._universes[k][j][i]
assert isinstance(u, Universe)
univs[u._id] = u
if self.outer is not None:
univs[self.outer._id] = self.outer
return univs
def get_all_nuclides(self):
"""Return all nuclides contained in the lattice
Returns
-------
nuclides : collections.OrderedDict
Dictionary whose keys are nuclide names and values are 2-tuples of
(nuclide, density)
"""
nuclides = OrderedDict()
# Get all unique Universes contained in each of the lattice cells
unique_universes = self.get_unique_universes()
# Append all Universes containing each cell to the dictionary
for universe_id, universe in unique_universes.items():
nuclides.update(universe.get_all_nuclides())
return nuclides
def get_all_cells(self):
"""Return all cells that are contained within the lattice
Returns
-------
cells : collections.OrderedDict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
"""
cells = OrderedDict()
unique_universes = self.get_unique_universes()
for universe_id, universe in unique_universes.items():
cells.update(universe.get_all_cells())
return cells
def get_all_materials(self):
"""Return all materials that are contained within the lattice
Returns
-------
materials : collections.OrderedDict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
"""
materials = OrderedDict()
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells()
for cell_id, cell in cells.items():
materials.update(cell.get_all_materials())
return materials
def get_all_universes(self):
"""Return all universes that are contained within the lattice
Returns
-------
universes : collections.OrderedDict
Dictionary whose keys are universe IDs and values are
:class:`Universe` instances
"""
# Initialize a dictionary of all Universes contained by the Lattice
# in each nested Universe level
all_universes = OrderedDict()
# Get all unique Universes contained in each of the lattice cells
unique_universes = self.get_unique_universes()
# Add the unique Universes filling each Lattice cell
all_universes.update(unique_universes)
# Append all Universes containing each cell to the dictionary
for universe_id, universe in unique_universes.items():
all_universes.update(universe.get_all_universes())
return all_universes
class RectLattice(Lattice):
"""A lattice consisting of rectangular prisms.
Parameters
----------
lattice_id : int, optional
Unique identifier for the lattice. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the lattice. If not specified, the name is the empty string.
Attributes
----------
id : int
Unique identifier for the lattice
name : str
Name of the lattice
dimension : Iterable of int
An array of two or three integers representing the number of lattice
cells in the x- and y- (and z-) directions, respectively.
lower_left : Iterable of float
The coordinates of the lower-left corner of the lattice. If the lattice
is two-dimensional, only the x- and y-coordinates are specified.
"""
def __init__(self, lattice_id=None, name=''):
super(RectLattice, self).__init__(lattice_id, name)
# Initialize Lattice class attributes
self._dimension = None
self._lower_left = None
self._offsets = None
def __eq__(self, other):
if not isinstance(other, RectLattice):
return False
elif not super(RectLattice, self).__eq__(other):
return False
elif self.dimension != other.dimension:
return False
elif self.lower_left != other.lower_left:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'RectLattice\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t',
self._dimension)
string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t',
self._lower_left)
string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch)
if self._outer is not None:
string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t',
self._outer._id)
else:
string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t',
self._outer)
string += '{0: <16}\n'.format('\tUniverses')
# Lattice nested Universe IDs - column major for Fortran
for i, universe in enumerate(np.ravel(self._universes)):
string += '{0} '.format(universe._id)
# Add a newline character every time we reach end of row of cells
if (i+1) % self._dimension[-1] == 0:
string += '\n'
string = string.rstrip('\n')
if self._offsets is not None:
string += '{0: <16}\n'.format('\tOffsets')
# Lattice cell offsets
for i, offset in enumerate(np.ravel(self._offsets)):
string += '{0} '.format(offset)
# Add a newline character when we reach end of row of cells
if (i+1) % self._dimension[-1] == 0:
string += '\n'
string = string.rstrip('\n')
return string
@property
def dimension(self):
return self._dimension
@property
def lower_left(self):
return self._lower_left
@property
def offsets(self):
return self._offsets
@dimension.setter
def dimension(self, dimension):
cv.check_type('lattice dimension', dimension, Iterable, Integral)
cv.check_length('lattice dimension', dimension, 2, 3)
for dim in dimension:
cv.check_greater_than('lattice dimension', dim, 0)
self._dimension = dimension
@lower_left.setter
def lower_left(self, lower_left):
cv.check_type('lattice lower left corner', lower_left, Iterable, Real)
cv.check_length('lattice lower left corner', lower_left, 2, 3)
self._lower_left = lower_left
@offsets.setter
def offsets(self, offsets):
cv.check_type('lattice offsets', offsets, Iterable)
self._offsets = offsets
@Lattice.pitch.setter
def pitch(self, pitch):
cv.check_type('lattice pitch', pitch, Iterable, Real)
cv.check_length('lattice pitch', pitch, 2, 3)
for dim in pitch:
cv.check_greater_than('lattice pitch', dim, 0.0)
self._pitch = pitch
def get_cell_instance(self, path, distribcell_index):
# Extract the lattice element from the path
next_index = path.index('-')
lat_id_indices = path[:next_index]
path = path[next_index+2:]
# Extract the lattice cell indices from the path
i1 = lat_id_indices.index('(')
i2 = lat_id_indices.index(')')
i = lat_id_indices[i1+1:i2]
lat_x = int(i.split(',')[0]) - 1
lat_y = int(i.split(',')[1]) - 1
lat_z = int(i.split(',')[2]) - 1
# For 2D Lattices
if len(self._dimension) == 2:
offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1]
offset += self._universes[lat_x][lat_y].get_cell_instance(path,
distribcell_index)
# For 3D Lattices
else:
offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1]
offset += self._universes[lat_z][lat_y][lat_x].get_cell_instance(
path, distribcell_index)
return offset
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
path = './lattice[@id=\'{0}\']'.format(self._id)
test = xml_element.find(path)
# If the element does contain the Lattice subelement, then return
if test is not None:
return
lattice_subelement = ET.Element("lattice")
lattice_subelement.set("id", str(self._id))
if len(self._name) > 0:
lattice_subelement.set("name", str(self._name))
# Export the Lattice cell pitch
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = ' '.join(map(str, self._pitch))
# Export the Lattice outer Universe (if specified)
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
self._outer.create_xml_subelement(xml_element)
# Export Lattice cell dimensions
dimension = ET.SubElement(lattice_subelement, "dimension")
dimension.text = ' '.join(map(str, self._dimension))
# Export Lattice lower left
lower_left = ET.SubElement(lattice_subelement, "lower_left")
lower_left.text = ' '.join(map(str, self._lower_left))
# Export the Lattice nested Universe IDs - column major for Fortran
universe_ids = '\n'
# 3D Lattices
if len(self._dimension) == 3:
for z in range(self._dimension[2]):
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[z][y][x]
# Append Universe ID to the Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
# Create XML subelement for this Universe
universe.create_xml_subelement(xml_element)
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# 2D Lattices
else:
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[y][x]
# Append Universe ID to Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
# Create XML subelement for this Universe
universe.create_xml_subelement(xml_element)
# Add newline character when we reach end of row of cells
universe_ids += '\n'
# Remove trailing newline character from Universe IDs string
universe_ids = universe_ids.rstrip('\n')
universes = ET.SubElement(lattice_subelement, "universes")
universes.text = universe_ids
# Append the XML subelement for this Lattice to the XML element
xml_element.append(lattice_subelement)
class HexLattice(Lattice):
"""A lattice consisting of hexagonal prisms.
Parameters
----------
lattice_id : int, optional
Unique identifier for the lattice. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the lattice. If not specified, the name is the empty string.
Attributes
----------
id : int
Unique identifier for the lattice
name : str
Name of the lattice
num_rings : int
Number of radial ring positions in the xy-plane
num_axial : int
Number of positions along the z-axis.
center : Iterable of float
Coordinates of the center of the lattice. If the lattice does not have
axial sections then only the x- and y-coordinates are specified
"""
def __init__(self, lattice_id=None, name=''):
super(HexLattice, self).__init__(lattice_id, name)
# Initialize Lattice class attributes
self._num_rings = None
self._num_axial = None
self._center = None
def __eq__(self, other):
if not isinstance(other, HexLattice):
return False
elif not super(HexLattice, self).__eq__(other):
return False
elif self.num_rings != other.num_rings:
return False
elif self.num_axial != other.num_axial:
return False
elif self.center != other.center:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'HexLattice\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings)
string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial)
string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t',
self._center)
string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch)
if self._outer is not None:
string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t',
self._outer._id)
else:
string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t',
self._outer)
string += '{0: <16}\n'.format('\tUniverses')
if self._num_axial is not None:
slices = [self._repr_axial_slice(x) for x in self._universes]
string += '\n'.join(slices)
else:
string += self._repr_axial_slice(self._universes)
return string
@property
def num_rings(self):
return self._num_rings
@property
def num_axial(self):
return self._num_axial
@property
def center(self):
return self._center
@num_rings.setter
def num_rings(self, num_rings):
cv.check_type('number of rings', num_rings, Integral)
cv.check_greater_than('number of rings', num_rings, 0)
self._num_rings = num_rings
@num_axial.setter
def num_axial(self, num_axial):
cv.check_type('number of axial', num_axial, Integral)
cv.check_greater_than('number of axial', num_axial, 0)
self._num_axial = num_axial
@center.setter
def center(self, center):
cv.check_type('lattice center', center, Iterable, Real)
cv.check_length('lattice center', center, 2, 3)
self._center = center
@Lattice.pitch.setter
def pitch(self, pitch):
cv.check_type('lattice pitch', pitch, Iterable, Real)
cv.check_length('lattice pitch', pitch, 1, 2)
for dim in pitch:
cv.check_greater_than('lattice pitch', dim, 0)
self._pitch = pitch
@Lattice.universes.setter
def universes(self, universes):
# Call Lattice.universes parent class setter property
Lattice.universes.fset(self, universes)
# NOTE: This routine assumes that the user creates a "ragged" list of
# lists, where each sub-list corresponds to one ring of Universes.
# The sub-lists are ordered from outermost ring to innermost ring.
# The Universes within each sub-list are ordered from the "top" in a
# clockwise fashion.
# Check to see if the given universes look like a 2D or a 3D array.
if isinstance(self._universes[0][0], Universe):
n_dims = 2
elif isinstance(self._universes[0][0][0], Universe):
n_dims = 3
else:
msg = 'HexLattice ID={0:d} does not appear to be either 2D or ' \
'3D. Make sure set_universes was given a two-deep or ' \
'three-deep iterable of universes.'.format(self._id)
raise RuntimeError(msg)
# Set the number of axial positions.
if n_dims == 3:
self.num_axial = len(self._universes)
else:
self._num_axial = None
# Set the number of rings and make sure this number is consistent for
# all axial positions.
if n_dims == 3:
self.num_rings = len(self._universes)
for rings in self._universes:
if len(rings) != self._num_rings:
msg = 'HexLattice ID={0:d} has an inconsistent number of ' \
'rings per axial positon'.format(self._id)
raise ValueError(msg)
else:
self.num_rings = len(self._universes)
# Make sure there are the correct number of elements in each ring.
if n_dims == 3:
for axial_slice in self._universes:
# Check the center ring.
if len(axial_slice[-1]) != 1:
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in the innermost ring. Only 1 element is ' \
'allowed in the innermost ring.'.format(self._id)
raise ValueError(msg)
# Check the outer rings.
for r in range(self._num_rings-1):
if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in ring number {1:d} (counting from the '\
'outermost ring). This ring should have {2:d} ' \
'elements.'.format(self._id, r,
6*(self._num_rings - 1 - r))
raise ValueError(msg)
else:
axial_slice = self._universes
# Check the center ring.
if len(axial_slice[-1]) != 1:
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in the innermost ring. Only 1 element is ' \
'allowed in the innermost ring.'.format(self._id)
raise ValueError(msg)
# Check the outer rings.
for r in range(self._num_rings-1):
if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
msg = 'HexLattice ID={0:d} has the wrong number of ' \
'elements in ring number {1:d} (counting from the '\
'outermost ring). This ring should have {2:d} ' \
'elements.'.format(self._id, r,
6*(self._num_rings - 1 - r))
raise ValueError(msg)
def create_xml_subelement(self, xml_element):
# Determine if XML element already contains subelement for this Lattice
path = './hex_lattice[@id=\'{0}\']'.format(self._id)
test = xml_element.find(path)
# If the element does contain the Lattice subelement, then return
if test is not None:
return
lattice_subelement = ET.Element("hex_lattice")
lattice_subelement.set("id", str(self._id))
if len(self._name) > 0:
lattice_subelement.set("name", str(self._name))
# Export the Lattice cell pitch
pitch = ET.SubElement(lattice_subelement, "pitch")
pitch.text = ' '.join(map(str, self._pitch))
# Export the Lattice outer Universe (if specified)
if self._outer is not None:
outer = ET.SubElement(lattice_subelement, "outer")
outer.text = '{0}'.format(self._outer._id)
self._outer.create_xml_subelement(xml_element)
lattice_subelement.set("n_rings", str(self._num_rings))
if self._num_axial is not None:
lattice_subelement.set("n_axial", str(self._num_axial))
# Export Lattice cell center
dimension = ET.SubElement(lattice_subelement, "center")
dimension.text = ' '.join(map(str, self._center))
# Export the Lattice nested Universe IDs.
# 3D Lattices
if self._num_axial is not None:
slices = []
for z in range(self._num_axial):
# Initialize the center universe.
universe = self._universes[z][-1][0]
universe.create_xml_subelement(xml_element)
# Initialize the remaining universes.
for r in range(self._num_rings-1):
for theta in range(6*(self._num_rings - 1 - r)):
universe = self._universes[z][r][theta]
universe.create_xml_subelement(xml_element)
# Get a string representation of the universe IDs.
slices.append(self._repr_axial_slice(self._universes[z]))
# Collapse the list of axial slices into a single string.
universe_ids = '\n'.join(slices)
# 2D Lattices
else:
# Initialize the center universe.
universe = self._universes[-1][0]
universe.create_xml_subelement(xml_element)
# Initialize the remaining universes.
for r in range(self._num_rings - 1):
for theta in range(6*(self._num_rings - 1 - r)):
universe = self._universes[r][theta]
universe.create_xml_subelement(xml_element)
# Get a string representation of the universe IDs.
universe_ids = self._repr_axial_slice(self._universes)
universes = ET.SubElement(lattice_subelement, "universes")
universes.text = '\n' + universe_ids
# Append the XML subelement for this Lattice to the XML element
xml_element.append(lattice_subelement)
def _repr_axial_slice(self, universes):
"""Return string representation for the given 2D group of universes.
The 'universes' argument should be a list of lists of universes where
each sub-list represents a single ring. The first list should be the
outer ring.
"""
# Find the largest universe ID and count the number of digits so we can
# properly pad the output string later.
largest_id = max([max([univ._id for univ in ring])
for ring in universes])
n_digits = len(str(largest_id))
pad = ' '*n_digits
id_form = '{: ^' + str(n_digits) + 'd}'
# Initialize the list for each row.
rows = [[] for i in range(1 + 4 * (self._num_rings-1))]
middle = 2 * (self._num_rings - 1)
# Start with the degenerate first ring.
universe = universes[-1][0]
rows[middle] = [id_form.format(universe._id)]
# Add universes one ring at a time.
for r in range(1, self._num_rings):
# r_prime increments down while r increments up.
r_prime = self._num_rings - 1 - r
theta = 0
y = middle + 2*r
# Climb down the top-right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 1
theta += 1
# Climb down the right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 2
theta += 1
# Climb down the bottom-right.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].append(id_form.format(universe._id))
# Translate the indices.
y -= 1
theta += 1
# Climb up the bottom-left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 1
theta += 1
# Climb up the left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 2
theta += 1
# Climb up the top-left.
for i in range(r):
# Add the universe.
universe = universes[r_prime][theta]
rows[y].insert(0, id_form.format(universe._id))
# Translate the indices.
y += 1
theta += 1
# Flip the rows and join each row into a single string.
rows = [pad.join(x) for x in rows[::-1]]
# Pad the beginning of the rows so they line up properly.
for y in range(self._num_rings - 1):
rows[y] = (self._num_rings - 1 - y)*pad + rows[y]
rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y]
for y in range(self._num_rings % 2, self._num_rings, 2):
rows[middle + y] = pad + rows[middle + y]
if y != 0:
rows[middle - y] = pad + rows[middle - y]
# Join the rows together and return the string.
universe_ids = '\n'.join(rows)
return universe_ids

21
openmc/material.py
View file

@ -270,7 +270,7 @@ class Material(object):
Parameters
----------
nuclide : str or openmc.nuclide.Nuclide
nuclide : str or openmc.Nuclide
Nuclide to add
percent : float
Atom or weight percent
@ -313,7 +313,7 @@ class Material(object):
Parameters
----------
nuclide : openmc.nuclide.Nuclide
nuclide : openmc.Nuclide
Nuclide to remove
"""
@ -332,7 +332,7 @@ class Material(object):
Parameters
----------
macroscopic : str or Macroscopic
macroscopic : str or openmc.Macroscopic
Macroscopic to add
"""
@ -371,7 +371,7 @@ class Material(object):
Parameters
----------
macroscopic : Macroscopic
macroscopic : openmc.Macroscopic
Macroscopic to remove
"""
@ -390,7 +390,7 @@ class Material(object):
Parameters
----------
element : openmc.element.Element
element : openmc.Element
Element to add
percent : float
Atom or weight percent
@ -429,7 +429,7 @@ class Material(object):
Parameters
----------
element : openmc.element.Element
element : openmc.Element
Element to remove
"""
@ -671,7 +671,7 @@ class MaterialsFile(object):
Parameters
----------
material : Material
material : openmc.Material
Material to add
"""
@ -688,7 +688,7 @@ class MaterialsFile(object):
Parameters
----------
materials : tuple or list of Material
materials : tuple or list of openmc.Material
Materials to add
"""
@ -706,7 +706,7 @@ class MaterialsFile(object):
Parameters
----------
material : Material
material : openmc.Material
Material to remove
"""
@ -723,9 +723,8 @@ class MaterialsFile(object):
material.make_isotropic_in_lab()
def _create_material_subelements(self):
subelement = ET.SubElement(self._materials_file, "default_xs")
if self._default_xs is not None:
subelement = ET.SubElement(self._materials_file, "default_xs")
subelement.text = self._default_xs
for material in self._materials:

16
openmc/mgxs/groups.py
View file

@ -24,7 +24,7 @@ class EnergyGroups(object):
----------
group_edges : Iterable of Real
The energy group boundaries [MeV]
num_groups : Integral
num_groups : int
The number of energy groups
"""
@ -86,7 +86,7 @@ class EnergyGroups(object):
Parameters
----------
energy : Real
energy : float
The energy of interest in MeV
Returns
@ -115,7 +115,7 @@ class EnergyGroups(object):
Parameters
----------
group : Integral
group : int
The energy group index, starting at 1 for the highest energies
Returns
@ -153,7 +153,7 @@ class EnergyGroups(object):
Returns
-------
ndarray
numpy.ndarray
The ndarray array indices for each energy group of interest
Raises
@ -200,7 +200,7 @@ class EnergyGroups(object):
Returns
-------
EnergyGroups
openmc.mgxs.EnergyGroups
A coarsened version of this EnergyGroups object.
Raises
@ -244,7 +244,7 @@ class EnergyGroups(object):
Parameters
----------
other : EnergyGroups
other : openmc.mgxs.EnergyGroups
EnergyGroups to compare with
Returns
@ -275,12 +275,12 @@ class EnergyGroups(object):
Parameters
----------
other : EnergyGroups
other : openmc.mgxs.EnergyGroups
EnergyGroups to merge with
Returns
-------
merged_groups : EnergyGroups
merged_groups : openmc.mgxs.EnergyGroups
EnergyGroups resulting from the merge
"""

22
openmc/mgxs/library.py
View file

@ -53,22 +53,22 @@ class Library(object):
The types of cross sections in the library (e.g., ['total', 'scatter'])
domain_type : {'material', 'cell', 'distribcell', 'universe'}
Domain type for spatial homogenization
domains : Iterable of Material, Cell or Universe
domains : Iterable of openmc.Material, openmc.Cell or openmc.Universe
The spatial domain(s) for which MGXS in the Library are computed
correction : 'P0' or None
correction : {'P0', None}
Apply the P0 correction to scattering matrices if set to 'P0'
energy_groups : EnergyGroups
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
tally_trigger : Trigger
tally_trigger : openmc.Trigger
An (optional) tally precision trigger given to each tally used to
compute the cross section
all_mgxs : OrderedDict
all_mgxs : collections.OrderedDict
MGXS objects keyed by domain ID and cross section type
sp_filename : str
The filename of the statepoint with tally data used to the
compute cross sections
keff : Real or None
The combined keff from the statepoint file with tally data used to
The combined keff from the statepoint file with tally data used to
compute cross sections (for eigenvalue calculations only)
name : str, optional
Name of the multi-group cross section library. Used as a label to
@ -308,7 +308,7 @@ class Library(object):
"""
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:
@ -350,7 +350,7 @@ class Library(object):
def add_to_tallies_file(self, tallies_file, merge=True):
"""Add all tallies from all MGXS objects to a tallies file.
NOTE: This assumes that build_library() has been called
NOTE: This assumes that :meth:`Library.build_library` has been called
Parameters
----------
@ -426,7 +426,7 @@ class Library(object):
----------
domain : Material or Cell or Universe or Integral
The material, cell, or universe object of interest (or its ID)
mgxs_type : {'total', 'transport', 'absorption', 'capture', 'fission', 'nu-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'chi'}
mgxs_type : {'total', 'transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'chi'}
The type of multi-group cross section object to return
Returns
@ -457,7 +457,7 @@ class Library(object):
break
else:
msg = 'Unable to find MGXS for {0} "{1}" in ' \
'library'.format(self.domain_type, domain)
'library'.format(self.domain_type, domain_id)
raise ValueError(msg)
else:
domain_id = domain.id
@ -537,7 +537,7 @@ class Library(object):
Returns
-------
Library
openmc.mgxs.Library
A new multi-group cross section library averaged across subdomains
Raises

107
openmc/mgxs/mgxs.py
View file

@ -59,18 +59,14 @@ class MGXS(object):
Parameters
----------
domain : Material or Cell or Universe
domain : openmc.Material or openmc.Cell or openmc.Universe
The domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe'}
The domain type for spatial homogenization
energy_groups : EnergyGroups
energy_groups : openmc.mgxs.EnergyGroups
The energy group structure for energy condensation
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
nuclides : Iterable of basestring
The user-specified nuclides to compute cross sections. If by_nuclide
is True but nuclides are not specified by the user, all nuclides in the
spatial domain will be used.
name : str, optional
Name of the multi-group cross section. Used as a label to identify
tallies in OpenMC 'tallies.xml' file.
@ -87,31 +83,33 @@ class MGXS(object):
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe'}
Domain type for spatial homogenization
energy_groups : EnergyGroups
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
tally_trigger : Trigger
tally_trigger : openmc.Trigger
An (optional) tally precision trigger given to each tally used to
compute the cross section
tallies : OrderedDict
tallies : collections.OrderedDict
OpenMC tallies needed to compute the multi-group cross section
rxn_rate_tally : Tally
rxn_rate_tally : openmc.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
xs_tally : openmc.Tally
Derived tally for the multi-group cross section. This attribute
is None unless the multi-group cross section has been computed.
num_subdomains : Integral
num_subdomains : int
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
tally data from a statepoint file).
num_nuclides : Integral
num_nuclides : int
The number of nuclides for which the multi-group cross section is
being tracked. This is unity if the by_nuclide attribute is False.
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.
nuclides : Iterable of str or 'sum'
The optional user-specified nuclides for which to compute cross
sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides
are not specified by the user, all nuclides in the spatial domain
are included. This attribute is 'sum' if by_nuclide is false.
sparse : bool
Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format
for compressed data storage
@ -336,11 +334,11 @@ class MGXS(object):
----------
mgxs_type : {'total', 'transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'chi'}
The type of multi-group cross section object to return
domain : Material or Cell or Universe
domain : openmc.Material or openmc.Cell or openmc.Universe
The domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe'}
The domain type for spatial homogenization
energy_groups : EnergyGroups
energy_groups : openmc.mgxs.EnergyGroups
The energy group structure for energy condensation
by_nuclide : bool
If true, computes cross sections for each nuclide in domain.
@ -351,7 +349,7 @@ class MGXS(object):
Returns
-------
MGXS
openmc.mgxs.MGXS
A subclass of the abstract MGXS class for the multi-group cross
section type requested by the user
@ -427,7 +425,7 @@ class MGXS(object):
Returns
-------
Real
float
The atomic number density (atom/b-cm) for the nuclide of interest
Raises
@ -466,7 +464,7 @@ class MGXS(object):
Returns
-------
ndarray of Real
numpy.ndarray of float
An array of the atomic number densities (atom/b-cm) for each of the
nuclides in the spatial domain
@ -514,11 +512,11 @@ class MGXS(object):
----------
scores : Iterable of str
Scores for each tally
all_filters : Iterable of tuple of Filter
all_filters : Iterable of tuple of openmc.Filter
Tuples of non-spatial domain filters for each tally
keys : Iterable of str
Key string used to store each tally in the tallies dictionary
estimator : {'analog' or 'tracklength'}
estimator : {'analog', 'tracklength'}
Type of estimator to use for each tally
"""
@ -537,27 +535,27 @@ class MGXS(object):
# Create each Tally needed to compute the multi group cross section
for score, key, filters in zip(scores, keys, all_filters):
self.tallies[key] = openmc.Tally(name=self.name)
self.tallies[key].add_score(score)
self.tallies[key].scores = [score]
self.tallies[key].estimator = estimator
self.tallies[key].add_filter(domain_filter)
self.tallies[key].filters = [domain_filter]
# If a tally trigger was specified, add it to each tally
if self.tally_trigger:
trigger_clone = copy.deepcopy(self.tally_trigger)
trigger_clone.add_score(score)
self.tallies[key].add_trigger(trigger_clone)
trigger_clone.scores = [score]
self.tallies[key].triggers.append(trigger_clone)
# Add all non-domain specific Filters (e.g., 'energy') to the Tally
for add_filter in filters:
self.tallies[key].add_filter(add_filter)
self.tallies[key].filters.append(add_filter)
# If this is a by-nuclide cross-section, add all nuclides to Tally
if self.by_nuclide and score != 'flux':
all_nuclides = self.domain.get_all_nuclides()
all_nuclides = self.get_all_nuclides()
for nuclide in all_nuclides:
self.tallies[key].add_nuclide(nuclide)
self.tallies[key].nuclides.append(nuclide)
else:
self.tallies[key].add_nuclide('total')
self.tallies[key].nuclides.append('total')
def _compute_xs(self):
"""Performs generic cleanup after a subclass' uses tally arithmetic to
@ -579,7 +577,7 @@ class MGXS(object):
self.xs_tally._nuclides = []
nuclides = self.get_all_nuclides()
for nuclide in nuclides:
self.xs_tally.add_nuclide(openmc.Nuclide(nuclide))
self.xs_tally.nuclides.append(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)
@ -686,7 +684,7 @@ class MGXS(object):
Returns
-------
ndarray
numpy.ndarray
A NumPy array of the multi-group cross section indexed in the order
each group, subdomain and nuclide is listed in the parameters.
@ -761,10 +759,9 @@ class MGXS(object):
# Reverse energies to align with increasing energy groups
xs = xs[:, ::-1, :]
# Eliminate trivial dimensions
xs = np.squeeze(xs)
xs = np.atleast_1d(xs)
# Eliminate trivial dimensions
xs = np.squeeze(xs)
xs = np.atleast_1d(xs)
return xs
def get_condensed_xs(self, coarse_groups):
@ -858,7 +855,7 @@ class MGXS(object):
Returns
-------
MGXS
openmc.mgxs.MGXS
A new MGXS averaged across the subdomains of interest
Raises
@ -910,13 +907,13 @@ class MGXS(object):
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
groups : list of Integral
groups : list of int
A list of energy group indices starting at 1 for the high energies
(e.g., [1, 2, 3]; default is [])
Returns
-------
MGXS
openmc.mgxs.MGXS
A new tally which encapsulates the subset of data requested for the
nuclide(s) and/or energy group(s) requested in the parameters.
@ -976,7 +973,7 @@ class MGXS(object):
Parameters
----------
other : MGXS
other : openmc.mgxs.MGXS
MGXS to check for merging
"""
@ -1013,12 +1010,12 @@ class MGXS(object):
Parameters
----------
other : MGXS
other : openmc.mgxs.MGXS
MGXS to merge with this one
Returns
-------
merged_mgxs : MGXS
merged_mgxs : openmc.mgxs.MGXS
Merged MGXS
"""
@ -1352,7 +1349,7 @@ class MGXS(object):
xs_type='macro', summary=None):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages the Tally.get_pandas_dataframe(...) method, but
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
renames the columns with terminology appropriate for cross section data.
Parameters
@ -1369,7 +1366,7 @@ class MGXS(object):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
summary : None or Summary
summary : None or openmc.Summary
An optional Summary object to be used to construct columns for
distribcell tally filters (default is None). The geometric
information in the Summary object is embedded into a multi-index
@ -1936,16 +1933,16 @@ class ScatterMatrixXS(MGXS):
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
in_groups : list of Integral
in_groups : list of int
A list of incoming energy group indices starting at 1 for the high
energies (e.g., [1, 2, 3]; default is [])
out_groups : list of Integral
out_groups : list of int
A list of outgoing energy group indices starting at 1 for the high
energies (e.g., [1, 2, 3]; default is [])
Returns
-------
MGXS
openmc.mgxs.MGXS
A new tally which encapsulates the subset of data requested for the
nuclide(s) and/or energy group(s) requested in the parameters.
@ -2313,7 +2310,7 @@ class Chi(MGXS):
super(Chi, self)._compute_xs()
# Add the coarse energy filter back to the nu-fission tally
nu_fission_in.add_filter(energy_filter)
nu_fission_in.filters.append(energy_filter)
return self._xs_tally
@ -2382,12 +2379,12 @@ class Chi(MGXS):
Parameters
----------
other : MGXS
other : openmc.mgxs.MGXS
MGXS to merge with this one
Returns
-------
merged_mgxs : MGXS
merged_mgxs : openmc.mgxs.MGXS
Merged MGXS
"""
@ -2455,7 +2452,7 @@ class Chi(MGXS):
Returns
-------
ndarray
numpy.ndarray
A NumPy array of the multi-group cross section indexed in the order
each group, subdomain and nuclide is listed in the parameters.
@ -2512,7 +2509,7 @@ class Chi(MGXS):
xs_tally = nu_fission_out / nu_fission_in
# Add the coarse energy filter back to the nu-fission tally
nu_fission_in.add_filter(energy_filter)
nu_fission_in.filters.append(energy_filter)
xs = xs_tally.get_values(filters=filters,
filter_bins=filter_bins, value=value)
@ -2563,7 +2560,7 @@ class Chi(MGXS):
xs_type='macro', summary=None):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages the Tally.get_pandas_dataframe(...) method, but
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
renames the columns with terminology appropriate for cross section data.
Parameters
@ -2580,7 +2577,7 @@ class Chi(MGXS):
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
summary : None or Summary
summary : None or openmc.Summary
An optional Summary object to be used to construct columns for
distribcell tally filters (default is None). The geometric
information in the Summary object is embedded into a multi-index

105
openmc/mgxs_library.py
View file

@ -24,7 +24,7 @@ def ndarray_to_string(arr):
Parameters
----------
arr : ndarray
arr : numpy.ndarray
Array to combine in to a string
Returns
@ -87,8 +87,9 @@ class XSdata(object):
----------
name : str, optional
Name of the mgxs data set.
representation : {'isotropic', 'angle'}
energy_groups : openmc.mgxs.EnergyGroups
Energygroup structure
representation : {'isotropic', 'angle'}, optional
Method used in generating the MGXS (isotropic or angle-dependent flux
weighting). Defaults to 'isotropic'
@ -99,10 +100,10 @@ class XSdata(object):
alias : str
Separate unique identifier for the xsdata object
kT : float
Temperature (in units of MeV) of this data set.
energy_groups : EnergyGroups
Temperature (in units of MeV).
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure
fissionable : boolean
fissionable : bool
Whether or not this is a fissionable data set.
scatt_type : {'legendre', 'histogram', or 'tabular'}
Angular distribution representation (legendre, histogram, or tabular)
@ -115,6 +116,85 @@ class XSdata(object):
Legendre polynomial form). Dict contains two keys: 'enable' and
'num_points'. 'enable' is a boolean and 'num_points' is the
number of points to use, if 'enable' is True.
num_azimuthal : int
Number of equal width angular bins that the azimuthal angular domain is
subdivided into. This only applies when ``representation`` is "angle".
num_polar : int
Number of equal width angular bins that the polar angular domain is
subdivided into. This only applies when ``representation`` is "angle".
total : numpy.ndarray
Group-wise total cross section ordered by increasing group index (i.e.,
fast to thermal). If ``representation`` is "isotropic", then the length
of this list should equal the number of groups described in the
``groups`` element. If ``representation`` is "angle", then the length
of this list should equal the number of groups times the number of
azimuthal angles times the number of polar angles, with the
inner-dimension being groups, intermediate-dimension being azimuthal
angles and outer-dimension being the polar angles.
absorption : numpy.ndarray
Group-wise absorption cross section ordered by increasing group index
(i.e., fast to thermal). If ``representation`` is "isotropic", then the
length of this list should equal the number of groups described in the
``groups`` attribute. If ``representation`` is "angle", then the length
of this list should equal the number of groups times the number of
azimuthal angles times the number of polar angles, with the
inner-dimension being groups, intermediate-dimension being azimuthal
angles and outer-dimension being the polar angles.
scatter : numpy.ndarray
Scattering moment matrices presented with the columns representing
incoming group and rows representing the outgoing group. That is,
down-scatter will be above the diagonal of the resultant matrix. This
matrix is repeated for every Legendre order (in order of increasing
orders) if ``scatt_type`` is "legendre"; otherwise, this matrix is
repeated for every bin of the histogram or tabular representation.
Finally, if ``representation`` is "angle", the above is repeated for
every azimuthal angle and every polar angle, in that order.
multiplicity : numpy.ndarray
Ratio of neutrons produced in scattering collisions to the neutrons
which undergo scattering collisions; that is, the multiplicity provides
the code with a scaling factor to account for neutrons being produced in
(n,xn) reactions. This information is assumed isotropic and therefore
does not need to be repeated for every Legendre moment or
histogram/tabular bin. This matrix follows the same arrangement as
described for the ``scatter`` attribute, with the exception of the data
needed to provide the scattering type information.
fission : numpy.ndarray
Group-wise fission cross section ordered by increasing group index
(i.e., fast to thermal). If ``representation`` is "isotropic", then the
length of this list should equal the number of groups described in the
``groups`` attribute. If ``representation`` is "angle", then the length
of this list should equal the number of groups times the number of
azimuthal angles times the number of polar angles, with the
inner-dimension being groups, intermediate-dimension being azimuthal
angles and outer-dimension being the polar angles.
k_fission : numpy.ndarray
Group-wise kappa-fission cross section ordered by increasing group index
(i.e., fast to thermal). If ``representation`` is "isotropic", then the
length of this list should equal the number of groups described in the
``groups`` attribute. If ``representation`` is "angle", then the length
of this list should equal the number of groups times the number of
azimuthal angles times the number of polar angles, with the
inner-dimension being groups, intermediate-dimension being azimuthal
angles and outer-dimension being the polar angles.
chi : numpy.ndarray
Group-wise fission spectra ordered by increasing group index (i.e., fast
to thermal). This attribute should be used if making the common
approximation that the fission spectra does not depend on incoming
energy. If the user does not wish to make this approximation, then this
should not be provided and this information included in the
``nu_fission`` element instead. If ``representation`` is "isotropic",
then the length of this list should equal the number of groups described
in the ``groups`` element. If ``representation`` is "angle", then the
length of this list should equal the number of groups times the number
of azimuthal angles times the number of polar angles, with the
inner-dimension being groups, intermediate-dimension being azimuthal
angles and outer-dimension being the polar angles.
nu_fission : numpy.ndarray
Group-wise fission production cross section vector (i.e., if ``chi`` is
provided), or is the group-wise fission production matrix. If providing
the vector, it should be ordered the same as the ``fission`` data. If
providing the matrix, it should be ordered the same as the
``multiplicity`` matrix.
"""
def __init__(self, name, energy_groups, representation="isotropic"):
@ -577,9 +657,7 @@ class MGXSLibraryFile(object):
Energy group structure.
inverse_velocities : Iterable of Real
Inverse of velocities, units of sec/cm
filename : str
XML file to write to.
xsdatas : Iterable of XSdata
xsdatas : Iterable of openmc.XSdata
Iterable of multi-Group cross section data objects
"""
@ -615,7 +693,7 @@ class MGXSLibraryFile(object):
Parameters
----------
xsdata : XSdata
xsdata : openmc.XSdata
MGXS information to add
"""
@ -638,7 +716,7 @@ class MGXSLibraryFile(object):
Parameters
----------
xsdatas : tuple or list of XSdata
xsdatas : tuple or list of openmc.XSdata
XSdatas to add
"""
@ -656,7 +734,7 @@ class MGXSLibraryFile(object):
Parameters
----------
xsdata : XSdata
xsdata : openmc.XSdata
XSdata to remove
"""
@ -717,6 +795,3 @@ class MGXSLibraryFile(object):
tree = ET.ElementTree(self._cross_sections_file)
tree.write(filename, xml_declaration=True,
encoding='utf-8', method="xml")

103
openmc/opencg_compatible.py
View file

@ -11,6 +11,7 @@ except ImportError:
import openmc
from openmc.region import Intersection
from openmc.surface import Halfspace
import openmc.checkvalue as cv
# A dictionary of all OpenMC Materials created
@ -79,10 +80,7 @@ def get_opencg_material(openmc_material):
"""
if not isinstance(openmc_material, openmc.Material):
msg = 'Unable to create an OpenCG Material from "{0}" ' \
'which is not an OpenMC Material'.format(openmc_material)
raise ValueError(msg)
cv.check_type('openmc_material', openmc_material, openmc.Material)
global OPENCG_MATERIALS
material_id = openmc_material.id
@ -119,10 +117,7 @@ def get_openmc_material(opencg_material):
"""
if not isinstance(opencg_material, opencg.Material):
msg = 'Unable to create an OpenMC Material from "{0}" ' \
'which is not an OpenCG Material'.format(opencg_material)
raise ValueError(msg)
cv.check_type('opencg_material', opencg_material, opencg.Material)
global OPENMC_MATERIALS
material_id = opencg_material.id
@ -165,10 +160,7 @@ def is_opencg_surface_compatible(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to check if OpenCG Surface is compatible' \
'since "{0}" is not a Surface'.format(opencg_surface)
raise ValueError(msg)
cv.check_type('opencg_surface', opencg_surface, opencg.Surface)
if opencg_surface.type in ['x-squareprism',
'y-squareprism', 'z-squareprism']:
@ -192,10 +184,7 @@ def get_opencg_surface(openmc_surface):
"""
if not isinstance(openmc_surface, openmc.Surface):
msg = 'Unable to create an OpenCG Surface from "{0}" ' \
'which is not an OpenMC Surface'.format(openmc_surface)
raise ValueError(msg)
cv.check_type('openmc_surface', openmc_surface, openmc.Surface)
global OPENCG_SURFACES
surface_id = openmc_surface.id
@ -278,10 +267,7 @@ def get_openmc_surface(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
cv.check_type('opencg_surface', opencg_surface, opencg.Surface)
global openmc_surface
surface_id = opencg_surface.id
@ -369,10 +355,7 @@ def get_compatible_opencg_surfaces(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
cv.check_type('opencg_surface', opencg_surface, opencg.Surface)
global OPENMC_SURFACES
surface_id = opencg_surface.id
@ -451,10 +434,7 @@ def get_opencg_cell(openmc_cell):
"""
if not isinstance(openmc_cell, openmc.Cell):
msg = 'Unable to create an OpenCG Cell from "{0}" which ' \
'is not an OpenMC Cell'.format(openmc_cell)
raise ValueError(msg)
cv.check_type('openmc_cell', openmc_cell, openmc.Cell)
global OPENCG_CELLS
cell_id = openmc_cell.id
@ -469,9 +449,9 @@ def get_opencg_cell(openmc_cell):
fill = openmc_cell.fill
if (openmc_cell.fill_type == 'material'):
if openmc_cell.fill_type == 'material':
opencg_cell.fill = get_opencg_material(fill)
elif (openmc_cell.fill_type == 'universe'):
elif openmc_cell.fill_type == 'universe':
opencg_cell.fill = get_opencg_universe(fill)
else:
opencg_cell.fill = get_opencg_lattice(fill)
@ -533,20 +513,10 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
OpenMC
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create compatible OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
elif not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create compatible OpenMC Cell since "{0}" is ' \
'not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
elif halfspace not in [-1, +1]:
msg = 'Unable to create compatible Cell since "{0}"' \
'is not a +/-1 halfspace'.format(halfspace)
raise ValueError(msg)
cv.check_type('opencg_cell', opencg_cell, opencg.Cell)
cv.check_type('opencg_surface', opencg_surface, opencg.Surface)
cv.check_value('halfspace', halfspace, (-1, +1))
# Initialize an empty list for the new compatible cells
compatible_cells = []
@ -575,7 +545,7 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
num_clones = 8
for clone_id in range(num_clones):
# Create a cloned OpenCG Cell with Surfaces compatible with OpenMC
# Create cloned OpenCG Cell with Surfaces compatible with OpenMC
clone = opencg_cell.clone()
compatible_cells.append(clone)
@ -641,10 +611,7 @@ def make_opencg_cells_compatible(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to make compatible OpenCG Cells for "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
cv.check_type('opencg_universe', opencg_universe, opencg.Universe)
# Check all OpenCG Cells in this Universe for compatibility with OpenMC
opencg_cells = opencg_universe.cells
@ -700,10 +667,7 @@ def get_openmc_cell(opencg_cell):
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create an OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
cv.check_type('opencg_cell', opencg_cell, opencg.Cell)
global OPENMC_CELLS
cell_id = opencg_cell.id
@ -718,9 +682,9 @@ def get_openmc_cell(opencg_cell):
fill = opencg_cell.fill
if (opencg_cell.type == 'universe'):
if opencg_cell.type == 'universe':
openmc_cell.fill = get_openmc_universe(fill)
elif (opencg_cell.type == 'lattice'):
elif opencg_cell.type == 'lattice':
openmc_cell.fill = get_openmc_lattice(fill)
else:
openmc_cell.fill = get_openmc_material(fill)
@ -764,10 +728,7 @@ def get_opencg_universe(openmc_universe):
"""
if not isinstance(openmc_universe, openmc.Universe):
msg = 'Unable to create an OpenCG Universe from "{0}" which ' \
'is not an OpenMC Universe'.format(openmc_universe)
raise ValueError(msg)
cv.check_type('openmc_universe', openmc_universe, openmc.Universe)
global OPENCG_UNIVERSES
universe_id = openmc_universe.id
@ -811,10 +772,7 @@ def get_openmc_universe(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to create an OpenMC Universe from "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
cv.check_type('opencg_universe', opencg_universe, opencg.Universe)
global OPENMC_UNIVERSES
universe_id = opencg_universe.id
@ -861,10 +819,7 @@ def get_opencg_lattice(openmc_lattice):
"""
if not isinstance(openmc_lattice, openmc.Lattice):
msg = 'Unable to create an OpenCG Lattice from "{0}" which ' \
'is not an OpenMC Lattice'.format(openmc_lattice)
raise ValueError(msg)
cv.check_type('openmc_lattice', openmc_lattice, openmc.Lattice)
global OPENCG_LATTICES
lattice_id = openmc_lattice.id
@ -958,10 +913,7 @@ def get_openmc_lattice(opencg_lattice):
"""
if not isinstance(opencg_lattice, opencg.Lattice):
msg = 'Unable to create an OpenMC Lattice from "{0}" which ' \
'is not an OpenCG Lattice'.format(opencg_lattice)
raise ValueError(msg)
cv.check_type('opencg_lattice', opencg_lattice, opencg.Lattice)
global OPENMC_LATTICES
lattice_id = opencg_lattice.id
@ -1032,10 +984,7 @@ def get_opencg_geometry(openmc_geometry):
"""
if not isinstance(openmc_geometry, openmc.Geometry):
msg = 'Unable to get OpenCG geometry from "{0}" which is ' \
'not an OpenMC Geometry object'.format(openmc_geometry)
raise ValueError(msg)
cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
# Clear dictionaries and auto-generated IDs
OPENMC_SURFACES.clear()
@ -1053,6 +1002,7 @@ def get_opencg_geometry(openmc_geometry):
opencg_geometry = opencg.Geometry()
opencg_geometry.root_universe = opencg_root_universe
opencg_geometry.initialize_cell_offsets()
opencg_geometry.assign_auto_ids()
return opencg_geometry
@ -1072,10 +1022,7 @@ def get_openmc_geometry(opencg_geometry):
"""
if not isinstance(opencg_geometry, opencg.Geometry):
msg = 'Unable to get OpenMC geometry from "{0}" which is ' \
'not an OpenCG Geometry object'.format(opencg_geometry)
raise ValueError(msg)
cv.check_type('opencg_geometry', opencg_geometry, opencg.Geometry)
# Deep copy the goemetry since it may be modified to make all Surfaces
# compatible with OpenMC's specifications

6
openmc/plots.py
View file

@ -275,7 +275,7 @@ class Plot(object):
The random number seed used to generate the color scheme
"""
cv.check_type('geometry', geometry, openmc.Geometry)
cv.check_type('seed', seed, Integral)
cv.check_greater_than('seed', seed, 1, equality=True)
@ -417,7 +417,7 @@ class PlotsFile(object):
Parameters
----------
plot : Plot
plot : openmc.Plot
Plot to add
"""
@ -433,7 +433,7 @@ class PlotsFile(object):
Parameters
----------
plot : Plot
plot : openmc.Plot
Plot to remove
"""

48
openmc/region.py
View file

@ -9,10 +9,11 @@ from openmc.checkvalue import check_type
class Region(object):
"""Region of space that can be assigned to a cell.
Region is an abstract base class that is inherited by Halfspace,
Intersection, Union, and Complement. Each of those respective classes are
typically not instantiated directly but rather are created through operators
of the Surface and Region classes.
Region is an abstract base class that is inherited by
:class:`openmc.Halfspace`, :class:`openmc.Intersection`,
:class:`openmc.Union`, and :class:`openmc.Complement`. Each of those
respective classes are typically not instantiated directly but rather are
created through operators of the Surface and Region classes.
"""
@ -201,11 +202,11 @@ class Intersection(Region):
"""Intersection of two or more regions.
Instances of Intersection are generally created via the __and__ operator
applied to two instances of Region. This is illustrated in the following
example:
applied to two instances of :class:`openmc.Region`. This is illustrated in
the following example:
>>> equator = openmc.surface.ZPlane(z0=0.0)
>>> earth = openmc.surface.Sphere(R=637.1e6)
>>> equator = openmc.ZPlane(z0=0.0)
>>> earth = openmc.Sphere(R=637.1e6)
>>> northern_hemisphere = -earth & +equator
>>> southern_hemisphere = -earth & -equator
>>> type(northern_hemisphere)
@ -213,12 +214,12 @@ class Intersection(Region):
Parameters
----------
*nodes
\*nodes
Regions to take the intersection of
Attributes
----------
nodes : tuple of Region
nodes : tuple of openmc.Region
Regions to take the intersection of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
@ -255,21 +256,22 @@ class Union(Region):
"""Union of two or more regions.
Instances of Union are generally created via the __or__ operator applied to
two instances of Region. This is illustrated in the following example:
two instances of :class:`openmc.Region`. This is illustrated in the
following example:
>>> s1 = openmc.surface.ZPlane(z0=0.0)
>>> s2 = openmc.surface.Sphere(R=637.1e6)
>>> s1 = openmc.ZPlane(z0=0.0)
>>> s2 = openmc.Sphere(R=637.1e6)
>>> type(-s2 | +s1)
<class 'openmc.region.Union'>
Parameters
----------
*nodes
\*nodes
Regions to take the union of
Attributes
----------
nodes : tuple of Region
nodes : tuple of openmc.Region
Regions to take the union of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
@ -305,13 +307,13 @@ class Union(Region):
class Complement(Region):
"""Complement of a region.
The Complement of an existing Region can be created by using the __invert__
operator as the following example demonstrates:
The Complement of an existing :class:`openmc.Region` can be created by using
the __invert__ operator as the following example demonstrates:
>>> xl = openmc.surface.XPlane(x0=-10.0)
>>> xr = openmc.surface.XPlane(x0=10.0)
>>> yl = openmc.surface.YPlane(y0=-10.0)
>>> yr = openmc.surface.YPlane(y0=10.0)
>>> xl = openmc.XPlane(x0=-10.0)
>>> xr = openmc.XPlane(x0=10.0)
>>> yl = openmc.YPlane(y0=-10.0)
>>> yr = openmc.YPlane(y0=10.0)
>>> inside_box = +xl & -xr & +yl & -yl
>>> outside_box = ~inside_box
>>> type(outside_box)
@ -319,12 +321,12 @@ class Complement(Region):
Parameters
----------
node : Region
node : openmc.Region
Region to take the complement of
Attributes
----------
node : Region
node : openmc.Region
Regions to take the complement of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box

136
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 import Nuclide
from openmc.source import Source
if sys.version_info[0] >= 3:
@ -37,7 +38,7 @@ 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 : Iterable of openmc.source.Source
source : Iterable of openmc.Source
Distribution of source sites in space, angle, and energy
output : dict
Dictionary indicating what files to output. Valid keys are 'summary',
@ -125,6 +126,8 @@ class SettingsFile(object):
Coordinates of the lower-left point of the UFS mesh
ufs_upper_right : tuple or list
Coordinates of the upper-right point of the UFS mesh
resonance_scattering : ResonanceScattering or iterable of ResonanceScattering
The elastic scattering model to use for resonant isotopes
"""
@ -205,6 +208,8 @@ class SettingsFile(object):
self._run_mode_subelement = None
self._source_element = None
self._resonance_scattering = None
@property
def run_mode(self):
return self._run_mode
@ -393,9 +398,13 @@ class SettingsFile(object):
def dd_count_interactions(self):
return self._dd_count_interactions
@property
def resonance_scattering(self):
return self._resonance_scattering
@run_mode.setter
def run_mode(self, run_mode):
if 'run_mode' not in ['eigenvalue', 'fixed source']:
if run_mode not in ['eigenvalue', 'fixed source']:
msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \
'and "fixed source" are supported."'.format(run_mode)
raise ValueError(msg)
@ -764,6 +773,16 @@ class SettingsFile(object):
self._dd_count_interactions = interactions
@resonance_scattering.setter
def resonance_scattering(self, res):
if isinstance(res, Iterable):
check_type('resonance_scattering', res, Iterable,
ResonanceScattering)
self._resonance_scattering = res
else:
check_type('resonance_scattering', res, ResonanceScattering)
self._resonance_scattering = [res]
def _create_run_mode_subelement(self):
if self.run_mode == 'eigenvalue':
@ -1043,6 +1062,17 @@ class SettingsFile(object):
subelement = ET.SubElement(element, "count_interactions")
subelement.text = str(self._dd_count_interactions).lower()
def _create_resonance_scattering_element(self):
if self.resonance_scattering is None: return
element = ET.SubElement(self._settings_file, "resonance_scattering")
for r in self.resonance_scattering:
if r.nuclide.name != r.nuclide_0K.name:
raise ValueError("The nuclide and nuclide_0K attributes of "
"a ResonantScattering object must have identical names.")
r.create_xml_subelement(element)
def export_to_xml(self):
"""Create a settings.xml file that can be used for a simulation.
@ -1079,6 +1109,7 @@ class SettingsFile(object):
self._create_track_subelement()
self._create_ufs_subelement()
self._create_dd_subelement()
self._create_resonance_scattering_element()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._settings_file)
@ -1087,3 +1118,104 @@ class SettingsFile(object):
tree = ET.ElementTree(self._settings_file)
tree.write("settings.xml", xml_declaration=True,
encoding='utf-8', method="xml")
class ResonanceScattering(object):
"""Specification of the elastic scattering model for resonant isotopes
Attributes
----------
nuclide : openmc.Nuclide
The nuclide affected by this resonance scattering treatment.
nuclide_0K : openmc.Nuclide
This should be the same isotope as the nuclide attribute above, but it
should have an xs attribute that identifies 0 Kelvin data.
method : str
The method used to sample outgoing scattering energies. Valid options
are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening
rejection correction), and 'WCM' (weight correction method).
E_min : float
The minimum energy above which the specified method is applied. By
default, CXS will be used below E_min.
E_max : float
The maximum energy below which the specified method is applied. By
default, the asymptotic target-at-rest model is applied above E_max.
"""
def __init__(self):
self._nuclide = None
self._nuclide_0K = None
self._method = None
self._E_min = None
self._E_max = None
@property
def nuclide(self):
return self._nuclide
@property
def nuclide_0K(self):
return self._nuclide_0K
@property
def method(self):
return self._method
@property
def E_min(self):
return self._E_min
@property
def E_max(self):
return self._E_max
@nuclide.setter
def nuclide(self, nuc):
check_type('nuclide', nuc, Nuclide)
if nuc.zaid == None: raise ValueError("The nuclide must have an "
"explicitly defined zaid attribute.")
self._nuclide = nuc
@nuclide_0K.setter
def nuclide_0K(self, nuc):
check_type('nuclide_0K', nuc, Nuclide)
if nuc.zaid == None: raise ValueError("The nuclide_0K must have an "
"explicitly defined zaid attribute.")
self._nuclide_0K = nuc
@method.setter
def method(self, m):
check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))
self._method = m
@E_min.setter
def E_min(self, E):
check_type('E_min', E, Real)
check_greater_than('E_min', E, 0, True)
self._E_min = E
@E_max.setter
def E_max(self, E):
check_type('E_max', E, Real)
check_greater_than('E_max', E, 0, True)
self._E_max = E
def create_xml_subelement(self, xml_element):
scatterer = ET.SubElement(xml_element, "scatterer")
subelement = ET.SubElement(scatterer, 'nuclide')
subelement.text = self.nuclide.name
if self.method is not None:
subelement = ET.SubElement(scatterer, 'method')
subelement.text = self.method
subelement = ET.SubElement(scatterer, 'xs_label')
subelement.text = str(self.nuclide.zaid) + '.' + str(self.nuclide.xs)
subelement = ET.SubElement(scatterer, 'xs_label_0K')
subelement.text = str(self.nuclide_0K.zaid) + '.' \
+ str(self.nuclide_0K.xs)
if self.E_min is not None:
subelement = ET.SubElement(scatterer, 'E_min')
subelement.text = str(self.E_min)
if self.E_max is not None:
subelement = ET.SubElement(scatterer, 'E_max')
subelement.text = str(self.E_max)

71
openmc/statepoint.py
View file

@ -18,59 +18,62 @@ class StatePoint(object):
----------
cmfd_on : bool
Indicate whether CMFD is active
cmfd_balance : ndarray
cmfd_balance : numpy.ndarray
Residual neutron balance for each batch
cmfd_dominance
Dominance ratio for each batch
cmfd_entropy : ndarray
cmfd_entropy : numpy.ndarray
Shannon entropy of CMFD fission source for each batch
cmfd_indices : ndarray
cmfd_indices : numpy.ndarray
Number of CMFD mesh cells and energy groups. The first three indices
correspond to the x-, y-, and z- spatial directions and the fourth index
is the number of energy groups.
cmfd_srccmp : ndarray
cmfd_srccmp : numpy.ndarray
Root-mean-square difference between OpenMC and CMFD fission source for
each batch
cmfd_src : ndarray
cmfd_src : numpy.ndarray
CMFD fission source distribution over all mesh cells and energy groups.
current_batch : Integral
current_batch : int
Number of batches simulated
date_and_time : str
Date and time when simulation began
entropy : ndarray
entropy : numpy.ndarray
Shannon entropy of fission source at each batch
gen_per_batch : Integral
Number of fission generations per batch
global_tallies : ndarray of compound datatype
global_tallies : numpy.ndarray of compound datatype
Global tallies for k-effective estimates and leakage. The compound
datatype has fields 'name', 'sum', 'sum_sq', 'mean', and 'std_dev'.
k_combined : list
Combined estimator for k-effective and its uncertainty
k_col_abs : Real
k_col_abs : float
Cross-product of collision and absorption estimates of k-effective
k_col_tra : Real
k_col_tra : float
Cross-product of collision and tracklength estimates of k-effective
k_abs_tra : Real
k_abs_tra : float
Cross-product of absorption and tracklength estimates of k-effective
k_generation : ndarray
k_generation : numpy.ndarray
Estimate of k-effective for each batch/generation
meshes : dict
Dictionary whose keys are mesh IDs and whose values are Mesh objects
n_batches : Integral
n_batches : int
Number of batches
n_inactive : Integral
n_inactive : int
Number of inactive batches
n_particles : Integral
n_particles : int
Number of particles per generation
n_realizations : Integral
n_realizations : int
Number of tally realizations
path : str
Working directory for simulation
run_mode : str
Simulation run mode, e.g. 'k-eigenvalue'
seed : Integral
runtime : dict
Dictionary whose keys are strings describing various runtime metrics
and whose values are time values in seconds.
seed : int
Pseudorandom number generator seed
source : ndarray of compound datatype
source : numpy.ndarray of compound datatype
Array of source sites. The compound datatype has fields 'wgt', 'xyz',
'uvw', and 'E' corresponding to the weight, position, direction, and
energy of the source site.
@ -88,7 +91,7 @@ class StatePoint(object):
TallyDerivative objects
version: tuple of Integral
Version of OpenMC
summary : None or openmc.summary.Summary
summary : None or openmc.Summary
A summary object if the statepoint has been linked with a summary file
"""
@ -104,8 +107,9 @@ class StatePoint(object):
raise IOError('{} is not a statepoint file.'.format(filename))
except AttributeError:
raise IOError('Could not read statepoint file. This most likely '
'means the statepoint file was produced by a different '
'version of OpenMC than the one you are using.')
'means the statepoint file was produced by a '
'different version of OpenMC than the one you are '
'using.')
if self._f['revision'].value != 15:
raise IOError('Statepoint file has a file revision of {} '
'which is not consistent with the revision this '
@ -315,6 +319,11 @@ class StatePoint(object):
def run_mode(self):
return self._f['run_mode'].value.decode()
@property
def runtime(self):
return {name: dataset.value
for name, dataset in self._f['runtime'].items()}
@property
def seed(self):
return self._f['seed'].value
@ -399,7 +408,7 @@ class StatePoint(object):
new_filter.mesh = self.meshes[key]
# Add Filter to the Tally
tally.add_filter(new_filter)
tally.filters.append(new_filter)
# Read Nuclide bins
nuclide_names = \
@ -408,7 +417,7 @@ class StatePoint(object):
# Add all Nuclides to the Tally
for name in nuclide_names:
nuclide = openmc.Nuclide(name.decode().strip())
tally.add_nuclide(nuclide)
tally.nuclides.append(nuclide)
scores = self._f['{0}{1}/score_bins'.format(
base, tally_key)].value
@ -435,7 +444,7 @@ class StatePoint(object):
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
tally.add_score(score)
tally.scores.append(score)
# Add Tally to the global dictionary of all Tallies
tally.sparse = self.sparse
@ -540,7 +549,7 @@ class StatePoint(object):
Returns
-------
tally : Tally
tally : openmc.Tally
A tally matching the specified criteria
Raises
@ -637,7 +646,7 @@ class StatePoint(object):
Parameters
----------
summary : Summary
summary : openmc.Summary
A Summary object.
Raises
@ -654,11 +663,13 @@ class StatePoint(object):
raise ValueError(msg)
for tally_id, tally in self.tallies.items():
# Get the Tally name from the summary file
tally.name = summary.tallies[tally_id].name
summary_tally = summary.tallies[tally_id]
tally.name = summary_tally.name
tally.with_summary = True
for tally_filter in tally.filters:
summary_filter = summary_tally.find_filter(tally_filter.type)
if tally_filter.type == 'surface':
surface_ids = []
for bin in tally_filter.bins:
@ -671,6 +682,10 @@ class StatePoint(object):
distribcell_ids.append(summary.cells[bin].id)
tally_filter.bins = distribcell_ids
if tally_filter.type == 'distribcell':
tally_filter.distribcell_paths = \
summary_filter.distribcell_paths
if tally_filter.type == 'universe':
universe_ids = []
for bin in tally_filter.bins:

40
openmc/stats/multivariate.py
View file

@ -22,12 +22,12 @@ class UnitSphere(object):
Parameters
----------
reference_uvw : Iterable of Real
reference_uvw : Iterable of float
Direction from which polar angle is measured
Attributes
----------
reference_uvw : Iterable of Real
reference_uvw : Iterable of float
Direction from which polar angle is measured
"""
@ -62,19 +62,19 @@ class PolarAzimuthal(UnitSphere):
Parameters
----------
mu : Univariate
mu : openmc.stats.Univariate
Distribution of the cosine of the polar angle
phi : Univariate
phi : openmc.stats.Univariate
Distribution of the azimuthal angle in radians
reference_uvw : Iterable of Real
reference_uvw : Iterable of float
Direction from which polar angle is measured. Defaults to the positive
z-direction.
Attributes
----------
mu : Univariate
mu : openmc.stats.Univariate
Distribution of the cosine of the polar angle
phi : Univariate
phi : openmc.stats.Univariate
Distribution of the azimuthal angle in radians
"""
@ -142,7 +142,7 @@ class Monodirectional(UnitSphere):
Parameters
----------
reference_uvw : Iterable of Real
reference_uvw : Iterable of float
Direction from which polar angle is measured. Defaults to the positive
x-direction.
@ -186,20 +186,20 @@ class CartesianIndependent(Spatial):
Parameters
----------
x : Univariate
x : openmc.stats.Univariate
Distribution of x-coordinates
y : Univariate
y : openmc.stats.Univariate
Distribution of y-coordinates
z : Univariate
z : openmc.stats.Univariate
Distribution of z-coordinates
Attributes
----------
x : Univariate
x : openmc.stats.Univariate
Distribution of x-coordinates
y : Univariate
y : openmc.stats.Univariate
Distribution of y-coordinates
z : Univariate
z : openmc.stats.Univariate
Distribution of z-coordinates
"""
@ -252,9 +252,9 @@ class Box(Spatial):
Parameters
----------
lower_left : Iterable of Real
lower_left : Iterable of float
Lower-left coordinates of cuboid
upper_right : Iterable of Real
upper_right : Iterable of float
Upper-right coordinates of cuboid
only_fissionable : bool, optional
Whether spatial sites should only be accepted if they occur in
@ -262,9 +262,9 @@ class Box(Spatial):
Attributes
----------
lower_left : Iterable of Real
lower_left : Iterable of float
Lower-left coordinates of cuboid
upper_right : Iterable of Real
upper_right : Iterable of float
Upper-right coordinates of cuboid
only_fissionable : bool, optional
Whether spatial sites should only be accepted if they occur in
@ -328,12 +328,12 @@ class Point(Spatial):
Parameters
----------
xyz : Iterable of Real
xyz : Iterable of float
Cartesian coordinates of location
Attributes
----------
xyz : Iterable of Real
xyz : Iterable of float
Cartesian coordinates of location
"""

16
openmc/stats/univariate.py
View file

@ -37,16 +37,16 @@ class Discrete(Univariate):
Parameters
----------
x : Iterable of Real
x : Iterable of float
Values of the random variable
p : Iterable of Real
p : Iterable of float
Discrete probability for each value
Attributes
----------
x : Iterable of Real
x : Iterable of float
Values of the random variable
p : Iterable of Real
p : Iterable of float
Discrete probability for each value
"""
@ -243,9 +243,9 @@ class Tabular(Univariate):
Parameters
----------
x : Iterable of Real
x : Iterable of float
Tabulated values of the random variable
p : Iterable of Real
p : Iterable of float
Tabulated probabilities
interpolation : {'histogram', 'linear-linear'}, optional
Indicate whether the density function is constant between tabulated
@ -253,9 +253,9 @@ class Tabular(Univariate):
Attributes
----------
x : Iterable of Real
x : Iterable of float
Tabulated values of the random variable
p : Iterable of Real
p : Iterable of float
Tabulated probabilities
interpolation : {'histogram', 'linear-linear'}, optional
Indicate whether the density function is constant between tabulated

20
openmc/summary.py
View file

@ -542,7 +542,7 @@ class Summary(object):
# 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)
tally.scores.append(score)
# Read filter metadata
num_filters = self._f['{0}/n_filters'.format(subbase)].value
@ -562,8 +562,14 @@ class Summary(object):
new_filter = openmc.Filter(filter_type, bins)
new_filter.num_bins = num_bins
# Read in distribcell paths
if filter_type == 'distribcell':
paths = self._f['{0}/paths'.format(subsubbase)][...]
paths = [str(path.decode()) for path in paths]
new_filter.distribcell_paths = paths
# Add Filter to the Tally
tally.add_filter(new_filter)
tally.filters.append(new_filter)
# Add Tally to the global dictionary of all Tallies
self.tallies[tally_id] = tally
@ -578,7 +584,7 @@ class Summary(object):
Returns
-------
material : openmc.material.Material
material : openmc.Material
Material with given id
"""
@ -599,7 +605,7 @@ class Summary(object):
Returns
-------
surface : openmc.surface.Surface
surface : openmc.Surface
Surface with given id
"""
@ -620,7 +626,7 @@ class Summary(object):
Returns
-------
cell : openmc.universe.Cell
cell : openmc.Cell
Cell with given id
"""
@ -641,7 +647,7 @@ class Summary(object):
Returns
-------
universe : openmc.universe.Universe
universe : openmc.Universe
Universe with given id
"""
@ -662,7 +668,7 @@ class Summary(object):
Returns
-------
lattice : openmc.universe.Lattice
lattice : openmc.Lattice
Lattice with given id
"""

49
openmc/surface.py
View file

@ -153,10 +153,10 @@ class Surface(object):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -278,8 +278,7 @@ class Plane(Surface):
class XPlane(Plane):
"""A plane perpendicular to the x axis, i.e. a surface of the form :math:`x -
x_0 = 0`
"""A plane perpendicular to the x axis of the form :math:`x - x_0 = 0`
Parameters
----------
@ -338,10 +337,10 @@ class XPlane(Plane):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -356,8 +355,7 @@ class XPlane(Plane):
class YPlane(Plane):
"""A plane perpendicular to the y axis, i.e. a surface of the form :math:`y -
y_0 = 0`
"""A plane perpendicular to the y axis of the form :math:`y - y_0 = 0`
Parameters
----------
@ -416,10 +414,10 @@ class YPlane(Plane):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -434,8 +432,7 @@ class YPlane(Plane):
class ZPlane(Plane):
"""A plane perpendicular to the z axis, i.e. a surface of the form :math:`z -
z_0 = 0`
"""A plane perpendicular to the z axis of the form :math:`z - z_0 = 0`
Parameters
----------
@ -494,10 +491,10 @@ class ZPlane(Plane):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -641,10 +638,10 @@ class XCylinder(Cylinder):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -740,10 +737,10 @@ class YCylinder(Cylinder):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -839,10 +836,10 @@ class ZCylinder(Cylinder):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -967,10 +964,10 @@ class Sphere(Surface):
Returns
-------
numpy.array
numpy.ndarray
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
numpy.ndarray
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
@ -1383,7 +1380,7 @@ class Halfspace(Region):
can be created from an existing Surface through the __neg__ and __pos__
operators, as the following example demonstrates:
>>> sphere = openmc.surface.Sphere(surface_id=1, R=10.0)
>>> sphere = openmc.Sphere(surface_id=1, R=10.0)
>>> inside_sphere = -sphere
>>> outside_sphere = +sphere
>>> type(inside_sphere)
@ -1391,18 +1388,18 @@ class Halfspace(Region):
Parameters
----------
surface : Surface
surface : openmc.Surface
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
Attributes
----------
surface : Surface
surface : openmc.Surface
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
bounding_box : tuple of numpy.array
bounding_box : tuple of numpy.ndarray
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""

401
openmc/tallies.py
View file

@ -1,14 +1,15 @@
from __future__ import division
from collections import Iterable, defaultdict
from collections import Iterable, MutableSequence, defaultdict
import copy
from functools import partial
import os
import pickle
import itertools
from numbers import Integral, Real
from xml.etree import ElementTree as ET
import sys
import warnings
from xml.etree import ElementTree as ET
import numpy as np
@ -18,10 +19,10 @@ from openmc.filter import _FILTER_TYPES
import openmc.checkvalue as cv
from openmc.clean_xml import *
if sys.version_info[0] >= 3:
basestring = str
# "Static" variable for auto-generated Tally IDs
AUTO_TALLY_ID = 10000
@ -32,6 +33,12 @@ AUTO_TALLY_ID = 10000
# specified axis.
_PRODUCT_TYPES = ['tensor', 'entrywise']
# The following indicate acceptable types when setting Tally.scores,
# Tally.nuclides, and Tally.filters
_SCORE_CLASSES = (basestring, CrossScore, AggregateScore)
_NUCLIDE_CLASSES = (basestring, Nuclide, CrossNuclide, AggregateNuclide)
_FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
def reset_auto_tally_id():
global AUTO_TALLY_ID
@ -44,7 +51,7 @@ class Tally(object):
Parameters
----------
tally_id : Integral, optional
tally_id : int, optional
Unique identifier for the tally. If none is specified, an identifier
will automatically be assigned
name : str, optional
@ -52,43 +59,43 @@ class Tally(object):
Attributes
----------
id : Integral
id : int
Unique identifier for the tally
name : str
Name of the tally
filters : list of openmc.filter.Filter
filters : list of openmc.Filter
List of specified filters for the tally
nuclides : list of openmc.nuclide.Nuclide
nuclides : list of openmc.Nuclide
List of nuclides to score results for
scores : list of str
List of defined scores, e.g. 'flux', 'fission', etc.
estimator : {'analog', 'tracklength', 'collision'}
Type of estimator for the tally
triggers : list of openmc.trigger.Trigger
triggers : list of openmc.Trigger
List of tally triggers
num_scores : Integral
num_scores : int
Total number of scores, accounting for the fact that a single
user-specified score, e.g. scatter-P3 or flux-Y2,2, might have multiple
bins
num_filter_bins : Integral
num_filter_bins : int
Total number of filter bins accounting for all filters
num_bins : Integral
num_bins : int
Total number of bins for the tally
shape : 3-tuple of Integral
The shape of the tally data array ordered as the number of filter bins,
shape : 3-tuple of int
The shape of the tally data array ordered as the number of filter bins,
nuclide bins and score bins
num_realizations : Integral
num_realizations : int
Total number of realizations
with_summary : bool
Whether or not a Summary has been linked
sum : ndarray
sum : numpy.ndarray
An array containing the sum of each independent realization for each bin
sum_sq : ndarray
sum_sq : numpy.ndarray
An array containing the sum of each independent realization squared for
each bin
mean : ndarray
mean : numpy.ndarray
An array containing the sample mean for each bin
std_dev : ndarray
std_dev : numpy.ndarray
An array containing the sample standard deviation for each bin
derived : bool
Whether or not the tally is derived from one or more other tallies
@ -104,12 +111,12 @@ class Tally(object):
# Initialize Tally class attributes
self.id = tally_id
self.name = name
self._filters = []
self._nuclides = []
self._scores = []
self._filters = cv.CheckedList(_FILTER_CLASSES, 'tally filters')
self._nuclides = cv.CheckedList(_NUCLIDE_CLASSES, 'tally nuclides')
self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores')
self._estimator = None
self._triggers = []
self._derivative = None
self._triggers = cv.CheckedList(Trigger, 'tally triggers')
self._num_realizations = 0
self._with_summary = False
@ -149,19 +156,19 @@ class Tally(object):
clone._filters = []
for self_filter in self.filters:
clone.add_filter(copy.deepcopy(self_filter, memo))
clone.filters.append(copy.deepcopy(self_filter, memo))
clone._nuclides = []
for nuclide in self.nuclides:
clone.add_nuclide(copy.deepcopy(nuclide, memo))
clone.nuclides.append(copy.deepcopy(nuclide, memo))
clone._scores = []
for score in self.scores:
clone.add_score(score)
clone.scores.append(score)
clone._triggers = []
for trigger in self.triggers:
clone.add_trigger(trigger)
clone.triggers.append(trigger)
memo[id(self)] = clone
@ -439,23 +446,30 @@ class Tally(object):
['analog', 'tracklength', 'collision'])
self._estimator = estimator
@triggers.setter
def triggers(self, triggers):
cv.check_type('tally triggers', triggers, MutableSequence)
self._triggers = cv.CheckedList(Trigger, 'tally triggers', triggers)
def add_trigger(self, trigger):
"""Add a tally trigger to the tally
.. deprecated:: 0.8
Use the Tally.triggers property directly, i.e.,
Tally.triggers.append(...)
Parameters
----------
trigger : openmc.trigger.Trigger
trigger : openmc.Trigger
Trigger to add
"""
if not isinstance(trigger, Trigger):
msg = 'Unable to add a tally trigger for Tally ID="{0}" to ' \
'since "{1}" is not a Trigger'.format(self.id, trigger)
raise ValueError(msg)
if trigger not in self.triggers:
self.triggers.append(trigger)
warnings.warn('Tally.add_trigger(...) has been deprecated and may be '
'removed in a future version. Tally triggers should be '
'defined using the triggers property directly.',
DeprecationWarning)
self.triggers.append(trigger)
@id.setter
def id(self, tally_id):
@ -482,9 +496,60 @@ class Tally(object):
cv.check_type('tally derivative', deriv, TallyDerivative)
self._derivative = deriv
@filters.setter
def filters(self, filters):
cv.check_type('tally filters', filters, MutableSequence)
# If the filter is already in the Tally, raise an error
for i, f in enumerate(filters[:-1]):
if f in filters[i+1:]:
msg = 'Unable to add a duplicate filter "{0}" to Tally ID="{1}" ' \
'since duplicate filters are not supported in the OpenMC ' \
'Python API'.format(f, self.id)
raise ValueError(msg)
self._filters = cv.CheckedList(_FILTER_CLASSES, 'tally filters', filters)
@nuclides.setter
def nuclides(self, nuclides):
cv.check_type('tally nuclides', nuclides, MutableSequence)
# If the nuclide is already in the Tally, raise an error
for i, nuclide in enumerate(nuclides[:-1]):
if nuclide in nuclides[i+1:]:
msg = 'Unable to add a duplicate nuclide "{0}" to Tally ID="{1}" ' \
'since duplicate nuclides are not supported in the OpenMC ' \
'Python API'.format(nuclide, self.id)
raise ValueError(msg)
self._nuclides = cv.CheckedList(_NUCLIDE_CLASSES, 'tally nuclides',
nuclides)
@scores.setter
def scores(self, scores):
cv.check_type('tally scores', scores, MutableSequence)
for i, score in enumerate(scores[:-1]):
# If the score is already in the Tally, raise an error
if score in scores[i+1:]:
msg = 'Unable to add a duplicate score "{0}" to Tally ID="{1}" ' \
'since duplicate scores are not supported in the OpenMC ' \
'Python API'.format(score, self.id)
raise ValueError(msg)
# If score is a string, strip whitespace
if isinstance(score, basestring):
scores[i] = score.strip()
self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores', scores)
def add_filter(self, new_filter):
"""Add a filter to the tally
.. deprecated:: 0.8
Use the Tally.filters property directly, i.e.,
Tally.filters.append(...)
Parameters
----------
new_filter : Filter, CrossFilter or AggregateFilter
@ -497,23 +562,19 @@ class Tally(object):
"""
if not isinstance(new_filter, (Filter, CrossFilter, AggregateFilter)):
msg = 'Unable to add Filter "{0}" to Tally ID="{1}" since it is ' \
'not a Filter object'.format(new_filter, self.id)
raise ValueError(msg)
# If the filter is already in the Tally, raise an error
if new_filter in self.filters:
msg = 'Unable to add a duplicate filter "{0}" to Tally ID="{1}" ' \
'since duplicate filters are not supported in the OpenMC ' \
'Python API'.format(new_filter, self.id)
raise ValueError(msg)
self._filters.append(new_filter)
warnings.warn('Tally.add_filter(...) has been deprecated and may be '
'removed in a future version. Tally filters should be '
'defined using the filters property directly.',
DeprecationWarning)
self.filters.append(new_filter)
def add_nuclide(self, nuclide):
"""Specify that scores for a particular nuclide should be accumulated
.. deprecated:: 0.8
Use the Tally.nuclides property directly, i.e.,
Tally.nuclides.append(...)
Parameters
----------
nuclide : str, Nuclide, CrossNuclide or AggregateNuclide
@ -526,24 +587,19 @@ class Tally(object):
"""
if not isinstance(nuclide, (basestring, Nuclide,
CrossNuclide, AggregateNuclide)):
msg = 'Unable to add nuclide "{0}" to Tally ID="{1}" since it is ' \
'not a Nuclide object'.format(nuclide)
raise ValueError(msg)
# If the nuclide is already in the Tally, raise an error
if nuclide in self.nuclides:
msg = 'Unable to add a duplicate nuclide "{0}" to Tally ID="{1}" ' \
'since duplicate nuclides are not supported in the OpenMC ' \
'Python API'.format(nuclide, self.id)
raise ValueError(msg)
self._nuclides.append(nuclide)
warnings.warn('Tally.add_nuclide(...) has been deprecated and may be '
'removed in a future version. Tally nuclides should be '
'defined using the nuclides property directly.',
DeprecationWarning)
self.nuclides.append(nuclide)
def add_score(self, score):
"""Specify a quantity to be scored
.. deprecated:: 0.8
Use the Tally.scores property directly, i.e.,
Tally.scores.append(...)
Parameters
----------
score : str, CrossScore or AggregateScore
@ -555,24 +611,11 @@ class Tally(object):
"""
if not isinstance(score, (basestring, CrossScore, AggregateScore)):
msg = 'Unable to add score "{0}" to Tally ID="{1}" since it is ' \
'not a string'.format(score, self.id)
raise ValueError(msg)
# If the score is already in the Tally, raise an error
if score in self.scores:
msg = 'Unable to add a duplicate score "{0}" to Tally ID="{1}" ' \
'since duplicate scores are not supported in the OpenMC ' \
'Python API'.format(score, self.id)
raise ValueError(msg)
# Normal score strings
if isinstance(score, basestring):
self._scores.append(score.strip())
# CrossScores and AggrgateScore
else:
self._scores.append(score)
warnings.warn('Tally.add_score(...) has been deprecated and may be '
'removed in a future version. Tally scores should be '
'defined using the scores property directly.',
DeprecationWarning)
self.scores.append(score)
@num_realizations.setter
def num_realizations(self, num_realizations):
@ -667,7 +710,7 @@ class Tally(object):
Parameters
----------
old_filter : openmc.filter.Filter
old_filter : openmc.Filter
Filter to remove
"""
@ -684,7 +727,7 @@ class Tally(object):
Parameters
----------
nuclide : openmc.nuclide.Nuclide
nuclide : openmc.Nuclide
Nuclide to remove
"""
@ -706,7 +749,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
Tally to check for mergeable filters
"""
@ -759,7 +802,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
Tally to check for mergeable nuclides
"""
@ -796,7 +839,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
Tally to check for mergeable scores
"""
@ -837,7 +880,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
Tally to check for merging
"""
@ -882,12 +925,12 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
Tally to merge with this one
Returns
-------
merged_tally : Tally
merged_tally : openmc.Tally
Merged tallies
"""
@ -939,7 +982,7 @@ class Tally(object):
# Add unique nuclides from other tally to merged tally
for nuclide in other.nuclides:
if nuclide not in merged_tally.nuclides:
merged_tally.add_nuclide(nuclide)
merged_tally.nuclides.append(nuclide)
# If two tallies can be merged along score bins
if merge_scores and not equal_scores:
@ -949,11 +992,11 @@ class Tally(object):
# Add unique scores from other tally to merged tally
for score in other.scores:
if score not in merged_tally.scores:
merged_tally.add_score(score)
merged_tally.scores.append(score)
# Add triggers from other tally to merged tally
for trigger in other.triggers:
merged_tally.add_trigger(trigger)
merged_tally.triggers.append(trigger)
# If results have not been read, then return tally for input generation
if self._results_read is None:
@ -1135,7 +1178,7 @@ class Tally(object):
Returns
-------
filter_found : openmc.filter.Filter
filter_found : openmc.Filter
Filter from this tally with matching type, or None if no matching
Filter is found
@ -1169,7 +1212,7 @@ class Tally(object):
----------
filter_type : str
The type of Filter (e.g., 'cell', 'energy', etc.)
filter_bin : Integral or tuple
filter_bin : int or tuple
The bin is an integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. The bin is an integer for the
cell instance ID for 'distribcell' Filters. The bin is a 2-tuple of
@ -1295,7 +1338,7 @@ class Tally(object):
Returns
-------
ndarray
numpy.ndarray
A NumPy array of the filter indices
"""
@ -1377,7 +1420,7 @@ class Tally(object):
Returns
-------
ndarray
numpy.ndarray
A NumPy array of the nuclide indices
"""
@ -1411,7 +1454,7 @@ class Tally(object):
Returns
-------
ndarray
numpy.ndarray
A NumPy array of the score indices
"""
@ -1473,7 +1516,7 @@ class Tally(object):
Returns
-------
float or ndarray
float or numpy.ndarray
A scalar or NumPy array of the Tally data indexed in the order
each filter, nuclide and score is listed in the parameters.
@ -1544,13 +1587,13 @@ class Tally(object):
Include columns with score bin information (default is True).
derivative : bool
Include columns with differential tally info (default is True).
summary : None or Summary
summary : None or openmc.Summary
An optional Summary object to be used to construct columns for
distribcell tally filters (default is None). The geometric
information in the Summary object is embedded into a Multi-index
column with a geometric "path" to each distribcell intance.
NOTE: This option requires the OpenCG Python package.
float_format : string
float_format : str
All floats in the DataFrame will be formatted using the given
format string before printing.
@ -1678,8 +1721,8 @@ class Tally(object):
The tally data in OpenMC is stored as a 3D array with the dimensions
corresponding to filters, nuclides and scores. As a result, tally data
can be opaque for a user to directly index (i.e., without use of the
Tally.get_values(...) method) since one must know how to properly use
can be opaque for a user to directly index (i.e., without use of
:meth:`openmc.Tally.get_values`) since one must know how to properly use
the number of bins and strides for each filter to index into the first
(filter) dimension.
@ -1699,7 +1742,7 @@ class Tally(object):
Returns
-------
ndarray
numpy.ndarray
The tally data array indexed by filters, nuclides and scores.
"""
@ -1877,7 +1920,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
The tally on the right hand side of the hybrid product
binary_op : {'+', '-', '*', '/', '^'}
The binary operation in the hybrid product
@ -1899,7 +1942,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new Tally that is the hybrid product with this one.
Raises
@ -2019,33 +2062,33 @@ class Tally(object):
# Add filters to the new tally
if filter_product == 'entrywise':
for self_filter in self_copy.filters:
new_tally.add_filter(self_filter)
new_tally.filters.append(self_filter)
else:
all_filters = [self_copy.filters, other_copy.filters]
for self_filter, other_filter in itertools.product(*all_filters):
new_filter = CrossFilter(self_filter, other_filter, binary_op)
new_tally.add_filter(new_filter)
new_tally.filters.append(new_filter)
# Add nuclides to the new tally
if nuclide_product == 'entrywise':
for self_nuclide in self_copy.nuclides:
new_tally.add_nuclide(self_nuclide)
new_tally.nuclides.append(self_nuclide)
else:
all_nuclides = [self_copy.nuclides, other_copy.nuclides]
for self_nuclide, other_nuclide in itertools.product(*all_nuclides):
new_nuclide = \
CrossNuclide(self_nuclide, other_nuclide, binary_op)
new_tally.add_nuclide(new_nuclide)
new_tally.nuclides.append(new_nuclide)
# Add scores to the new tally
if score_product == 'entrywise':
for self_score in self_copy.scores:
new_tally.add_score(self_score)
new_tally.scores.append(self_score)
else:
all_scores = [self_copy.scores, other_copy.scores]
for self_score, other_score in itertools.product(*all_scores):
new_score = CrossScore(self_score, other_score, binary_op)
new_tally.add_score(new_score)
new_tally.scores.append(new_score)
# Update the new tally's filter strides
new_tally._update_filter_strides()
@ -2077,7 +2120,7 @@ class Tally(object):
Parameters
----------
other : Tally
other : openmc.Tally
The tally to outer product with this tally
filter_product : {'entrywise'}
The type of product to be performed between filter data. Currently,
@ -2108,14 +2151,14 @@ class Tally(object):
filter_copy = copy.deepcopy(other_filter)
other._mean = np.repeat(other.mean, filter_copy.num_bins, axis=0)
other._std_dev = np.repeat(other.std_dev, filter_copy.num_bins, axis=0)
other.add_filter(filter_copy)
other.filters.append(filter_copy)
# Add filters present in other but not in self to self
for self_filter in self_missing_filters:
filter_copy = copy.deepcopy(self_filter)
self._mean = np.repeat(self.mean, filter_copy.num_bins, axis=0)
self._std_dev = np.repeat(self.std_dev, filter_copy.num_bins, axis=0)
self.add_filter(filter_copy)
self.filters.append(filter_copy)
# Align other filters with self filters
for i, self_filter in enumerate(self.filters):
@ -2138,7 +2181,7 @@ class Tally(object):
np.tile(other.std_dev, (1, self.num_nuclides, 1))
# Add nuclides to each tally such that each tally contains the complete
# set of nuclides necessary to perform an entrywise product. New
# set of nuclides necessary to perform an entrywise product. New
# nuclides added to a tally will have all their scores set to zero.
else:
@ -2154,7 +2197,7 @@ class Tally(object):
np.insert(other.mean, other.num_nuclides, 0, axis=1)
other._std_dev = \
np.insert(other.std_dev, other.num_nuclides, 0, axis=1)
other.add_nuclide(nuclide)
other.nuclides.append(nuclide)
# Add nuclides present in other but not in self to self
for nuclide in self_missing_nuclides:
@ -2162,7 +2205,7 @@ class Tally(object):
np.insert(self.mean, self.num_nuclides, 0, axis=1)
self._std_dev = \
np.insert(self.std_dev, self.num_nuclides, 0, axis=1)
self.add_nuclide(nuclide)
self.nuclides.append(nuclide)
# Align other nuclides with self nuclides
for i, nuclide in enumerate(self.nuclides):
@ -2195,13 +2238,13 @@ class Tally(object):
for score in other_missing_scores:
other._mean = np.insert(other.mean, other.num_scores, 0, axis=2)
other._std_dev = np.insert(other.std_dev, other.num_scores, 0, axis=2)
other.add_score(score)
other.scores.append(score)
# Add scores present in other but not in self to self
for score in self_missing_scores:
self._mean = np.insert(self.mean, self.num_scores, 0, axis=2)
self._std_dev = np.insert(self.std_dev, self.num_scores, 0, axis=2)
self.add_score(score)
self.scores.append(score)
# Align other scores with self scores
for i, score in enumerate(self.scores):
@ -2459,12 +2502,12 @@ class Tally(object):
Parameters
----------
other : Tally or Real
other : openmc.Tally or float
The tally or scalar value to add to this tally
Returns
-------
Tally
openmc.Tally
A new derived tally which is the sum of this tally and the other
tally or scalar value in the addition.
@ -2499,12 +2542,9 @@ class Tally(object):
new_tally.with_summary = self.with_summary
new_tally.num_realization = self.num_realizations
for self_filter in self.filters:
new_tally.add_filter(self_filter)
for nuclide in self.nuclides:
new_tally.add_nuclide(nuclide)
for score in self.scores:
new_tally.add_score(score)
new_tally.filters = copy.deepcopy(self.filters)
new_tally.nuclides = copy.deepcopy(self.nuclides)
new_tally.scores = copy.deepcopy(self.scores)
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
@ -2534,12 +2574,12 @@ class Tally(object):
Parameters
----------
other : Tally or Real
other : openmc.Tally or float
The tally or scalar value to subtract from this tally
Returns
-------
Tally
openmc.Tally
A new derived tally which is the difference of this tally and the
other tally or scalar value in the subtraction.
@ -2573,12 +2613,9 @@ class Tally(object):
new_tally.with_summary = self.with_summary
new_tally.num_realization = self.num_realizations
for self_filter in self.filters:
new_tally.add_filter(self_filter)
for nuclide in self.nuclides:
new_tally.add_nuclide(nuclide)
for score in self.scores:
new_tally.add_score(score)
new_tally.filters = copy.deepcopy(self.filters)
new_tally.nuclides = copy.deepcopy(self.nuclides)
new_tally.scores = copy.deepcopy(self.scores)
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
@ -2609,12 +2646,12 @@ class Tally(object):
Parameters
----------
other : Tally or Real
other : openmc.Tally or float
The tally or scalar value to multiply with this tally
Returns
-------
Tally
openmc.Tally
A new derived tally which is the product of this tally and the
other tally or scalar value in the multiplication.
@ -2648,12 +2685,9 @@ class Tally(object):
new_tally.with_summary = self.with_summary
new_tally.num_realization = self.num_realizations
for self_filter in self.filters:
new_tally.add_filter(self_filter)
for nuclide in self.nuclides:
new_tally.add_nuclide(nuclide)
for score in self.scores:
new_tally.add_score(score)
new_tally.filters = copy.deepcopy(self.filters)
new_tally.nuclides = copy.deepcopy(self.nuclides)
new_tally.scores = copy.deepcopy(self.scores)
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
@ -2684,12 +2718,12 @@ class Tally(object):
Parameters
----------
other : Tally or Real
other : openmc.Tally or float
The tally or scalar value to divide this tally by
Returns
-------
Tally
openmc.Tally
A new derived tally which is the dividend of this tally and the
other tally or scalar value in the division.
@ -2723,12 +2757,9 @@ class Tally(object):
new_tally.with_summary = self.with_summary
new_tally.num_realization = self.num_realizations
for self_filter in self.filters:
new_tally.add_filter(self_filter)
for nuclide in self.nuclides:
new_tally.add_nuclide(nuclide)
for score in self.scores:
new_tally.add_score(score)
new_tally.filters = copy.deepcopy(self.filters)
new_tally.nuclides = copy.deepcopy(self.nuclides)
new_tally.scores = copy.deepcopy(self.scores)
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
@ -2762,12 +2793,12 @@ class Tally(object):
Parameters
----------
power : Tally or Real
power : openmc.Tally or float
The tally or scalar value exponent
Returns
-------
Tally
openmc.Tally
A new derived tally which is this tally raised to the power of the
other tally or scalar value in the exponentiation.
@ -2802,12 +2833,9 @@ class Tally(object):
new_tally.with_summary = self.with_summary
new_tally.num_realization = self.num_realizations
for self_filter in self.filters:
new_tally.add_filter(self_filter)
for nuclide in self.nuclides:
new_tally.add_nuclide(nuclide)
for score in self.scores:
new_tally.add_score(score)
new_tally.filters = copy.deepcopy(self.filters)
new_tally.nuclides = copy.deepcopy(self.nuclides)
new_tally.scores = copy.deepcopy(self.scores)
# If original tally was sparse, sparsify the exponentiated tally
new_tally.sparse = self.sparse
@ -2826,12 +2854,12 @@ class Tally(object):
Parameters
----------
other : Integer or Real
other : float
The scalar value to add to this tally
Returns
-------
Tally
openmc.Tally
A new derived tally of this tally added with the scalar value.
"""
@ -2845,12 +2873,12 @@ class Tally(object):
Parameters
----------
other : Integer or Real
other : float
The scalar value to subtract this tally from
Returns
-------
Tally
openmc.Tally
A new derived tally of this tally subtracted from the scalar value.
"""
@ -2864,12 +2892,12 @@ class Tally(object):
Parameters
----------
other : Integer or Real
other : float
The scalar value to multiply with this tally
Returns
-------
Tally
openmc.Tally
A new derived tally of this tally multiplied by the scalar value.
"""
@ -2883,12 +2911,12 @@ class Tally(object):
Parameters
----------
other : Integer or Real
other : float
The scalar value to divide by this tally
Returns
-------
Tally
openmc.Tally
A new derived tally of the scalar value divided by this tally.
"""
@ -2900,7 +2928,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new derived tally which is the absolute value of this tally.
"""
@ -2914,7 +2942,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new derived tally which is the negated value of this tally.
"""
@ -2956,7 +2984,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new tally which encapsulates the subset of data requested in the
order each filter, nuclide and score is listed in the parameters.
@ -3079,7 +3107,7 @@ class Tally(object):
filter_type : str
A filter type string (e.g., 'cell', 'energy') corresponding to the
filter bins to sum across
filter_bins : Iterable of Integral or tuple
filter_bins : Iterable of int or tuple
A list of the filter bins corresponding to the filter_type parameter
Each bin in the list is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
@ -3097,7 +3125,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new tally which encapsulates the sum of data requested.
"""
@ -3148,11 +3176,11 @@ class Tally(object):
if not remove_filter:
filter_sum = \
AggregateFilter(self_filter, [tuple(filter_bins)], 'sum')
tally_sum.add_filter(filter_sum)
tally_sum.filters.append(filter_sum)
# Add a copy of each filter not summed across to the tally sum
else:
tally_sum.add_filter(copy.deepcopy(self_filter))
tally_sum.filters.append(copy.deepcopy(self_filter))
# Add a copy of this tally's filters to the tally sum
else:
@ -3170,7 +3198,7 @@ class Tally(object):
# Add AggregateNuclide to the tally sum
nuclide_sum = AggregateNuclide(nuclides, 'sum')
tally_sum.add_nuclide(nuclide_sum)
tally_sum.nuclides.append(nuclide_sum)
# Add a copy of this tally's nuclides to the tally sum
else:
@ -3188,7 +3216,7 @@ class Tally(object):
# Add AggregateScore to the tally sum
score_sum = AggregateScore(scores, 'sum')
tally_sum.add_score(score_sum)
tally_sum.scores.append(score_sum)
# Add a copy of this tally's scores to the tally sum
else:
@ -3227,7 +3255,7 @@ class Tally(object):
filter_type : str
A filter type string (e.g., 'cell', 'energy') corresponding to the
filter bins to average across
filter_bins : Iterable of Integral or tuple
filter_bins : Iterable of int or tuple
A list of the filter bins corresponding to the filter_type parameter
Each bin in the list is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
@ -3245,7 +3273,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new tally which encapsulates the average of data requested.
"""
@ -3297,11 +3325,11 @@ class Tally(object):
if not remove_filter:
filter_sum = \
AggregateFilter(self_filter, [tuple(filter_bins)], 'avg')
tally_avg.add_filter(filter_sum)
tally_avg.filters.append(filter_sum)
# Add a copy of each filter not averaged across to the tally avg
else:
tally_avg.add_filter(copy.deepcopy(self_filter))
tally_avg.filters.append(copy.deepcopy(self_filter))
# Add a copy of this tally's filters to the tally avg
else:
@ -3320,7 +3348,7 @@ class Tally(object):
# Add AggregateNuclide to the tally avg
nuclide_avg = AggregateNuclide(nuclides, 'avg')
tally_avg.add_nuclide(nuclide_avg)
tally_avg.nuclides.append(nuclide_avg)
# Add a copy of this tally's nuclides to the tally avg
else:
@ -3339,7 +3367,7 @@ class Tally(object):
# Add AggregateScore to the tally avg
score_sum = AggregateScore(scores, 'avg')
tally_avg.add_score(score_sum)
tally_avg.scores.append(score_sum)
# Add a copy of this tally's scores to the tally avg
else:
@ -3378,7 +3406,7 @@ class Tally(object):
Returns
-------
Tally
openmc.Tally
A new derived Tally with data diagaonalized along the new filter.
"""
@ -3392,7 +3420,7 @@ class Tally(object):
# Add the new filter to a copy of this Tally
new_tally = copy.deepcopy(self)
new_tally.add_filter(new_filter)
new_tally.filters.append(new_filter)
# Determine "base" indices along the new "diagonal", and the factor
# by which the "base" indices should be repeated to account for all
@ -3454,9 +3482,8 @@ class TalliesFile(object):
Parameters
----------
tally : Tally
tally : openmc.Tally
Tally to add to file
merge : bool
Indicate whether the tally should be merged with an existing tally,
if possible. Defaults to False.
@ -3493,7 +3520,7 @@ class TalliesFile(object):
Parameters
----------
tally : Tally
tally : openmc.Tally
Tally to remove
"""
@ -3529,7 +3556,7 @@ class TalliesFile(object):
Parameters
----------
mesh : openmc.mesh.Mesh
mesh : openmc.Mesh
Mesh to add to the file
"""
@ -3545,7 +3572,7 @@ class TalliesFile(object):
Parameters
----------
mesh : openmc.mesh.Mesh
mesh : openmc.Mesh
Mesh to remove from the file
"""

38
openmc/trigger.py
View file

@ -1,8 +1,10 @@
from numbers import Real
from xml.etree import ElementTree as ET
import sys
import warnings
from collections import Iterable
from openmc.checkvalue import check_type, check_value
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
@ -46,9 +48,7 @@ class Trigger(object):
clone._trigger_type = self._trigger_type
clone._threshold = self._threshold
clone._scores = []
for score in self._scores:
clone.add_score(score)
clone.scores = self.scores
memo[id(self)] = clone
@ -88,15 +88,26 @@ class Trigger(object):
@trigger_type.setter
def trigger_type(self, trigger_type):
check_value('tally trigger type', trigger_type,
cv.check_value('tally trigger type', trigger_type,
['variance', 'std_dev', 'rel_err'])
self._trigger_type = trigger_type
@threshold.setter
def threshold(self, threshold):
check_type('tally trigger threshold', threshold, Real)
cv.check_type('tally trigger threshold', threshold, Real)
self._threshold = threshold
@scores.setter
def scores(self, scores):
cv.check_type('trigger scores', scores, Iterable, basestring)
# Set scores making sure not to have duplicates
self._scores = []
for score in scores:
if score not in self._scores:
self._scores.append(score)
def add_score(self, score):
"""Add a score to the list of scores to be checked against the trigger.
@ -107,16 +118,11 @@ class Trigger(object):
"""
if not isinstance(score, basestring):
msg = 'Unable to add score "{0}" to tally trigger since ' \
'it is not a string'.format(score)
raise ValueError(msg)
# If the score is already in the Tally, don't add it again
if score in self._scores:
return
else:
self._scores.append(score)
warnings.warn('Trigger.add_score(...) has been deprecated and may be '
'removed in a future version. Tally trigger scores should '
'be defined using the scores property directly.',
DeprecationWarning)
self.scores.append(score)
def get_trigger_xml(self, element):
"""Return XML representation of the trigger

1334
openmc/universe.py
View file

File diff suppressed because it is too large Load diff

2
setup.py
View file

@ -36,7 +36,7 @@ if have_setuptools:
# Optional dependencies
'extras_require': {
'pandas': ['pandas'],
'pandas': ['pandas>=0.17.0'],
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],
'validate': ['lxml']

745
src/ace.F90
View file

File diff suppressed because it is too large Load diff

58
src/ace_header.F90
View file

@ -1,58 +0,0 @@
module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use secondary_header, only: SecondaryDistribution, AngleEnergyContainer
use stl_vector, only: VectorInt
implicit none
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
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
real(8), allocatable :: sigma(:) ! Cross section values
type(SecondaryDistribution) :: secondary
contains
procedure :: clear => reaction_clear ! Deallocates Reaction
end type Reaction
!===============================================================================
! URRDATA contains probability tables for the unresolved resonance range.
!===============================================================================
type UrrData
integer :: n_energy ! # of incident neutron energies
integer :: n_prob ! # of probabilities
integer :: interp ! inteprolation (2=lin-lin, 5=log-log)
integer :: inelastic_flag ! inelastic competition flag
integer :: absorption_flag ! other absorption flag
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
end type UrrData
contains
!===============================================================================
! REACTION_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine reaction_clear(this)
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
end subroutine reaction_clear
end module ace_header

29
src/angleenergy_header.F90 Normal file
View file

@ -0,0 +1,29 @@
module angleenergy_header
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
! distribution that is a function of incoming energy. Each derived type must
! implement a sample() subroutine that returns an outgoing energy and scattering
! cosine given an incoming energy.
!===============================================================================
type, abstract :: AngleEnergy
contains
procedure(angleenergy_sample_), deferred :: sample
end type AngleEnergy
abstract interface
subroutine angleenergy_sample_(this, E_in, E_out, mu)
import AngleEnergy
class(AngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in
real(8), intent(out) :: E_out
real(8), intent(out) :: mu
end subroutine angleenergy_sample_
end interface
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
end module angleenergy_header

2
src/cmfd_input.F90
View file

@ -70,7 +70,7 @@ contains
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
! CMFD is optional unless it is in on from settings
if (cmfd_on) then
if (cmfd_run) then
call fatal_error("No CMFD XML file, '" // trim(filename) // "' does not&
& exist!")
end if

15
src/constants.F90
View file

@ -223,6 +223,12 @@ module constants
NU_POLYNOMIAL = 1, & ! Nu values given by polynomial
NU_TABULAR = 2 ! Nu values given by tabular distribution
! Secondary particle emission type
integer, parameter :: &
EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle
EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle
EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed)
! Cross section filetypes
integer, parameter :: &
ASCII = 1, & ! ASCII cross section file
@ -370,10 +376,11 @@ module constants
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS
integer, parameter :: N_STREAMS = 3
integer, parameter :: STREAM_TRACKING = 1
integer, parameter :: STREAM_TALLIES = 2
integer, parameter :: STREAM_SOURCE = 3
integer, parameter :: N_STREAMS = 4
integer, parameter :: STREAM_TRACKING = 1
integer, parameter :: STREAM_TALLIES = 2
integer, parameter :: STREAM_SOURCE = 3
integer, parameter :: STREAM_URR_PTABLE = 4
! ============================================================================
! MISCELLANEOUS CONSTANTS

236
src/cross_section.F90
View file

@ -1,16 +1,14 @@
module cross_section
use ace_header, only: Reaction, UrrData
use constants
use energy_grid, only: grid_method, log_spacing
use error, only: fatal_error
use fission, only: nu_total
use global
use list_header, only: ListElemInt
use material_header, only: Material
use nuclide_header
use particle_header, only: Particle
use random_lcg, only: prn
use random_lcg, only: prn, future_prn, prn_set_stream
use sab_header, only: SAlphaBeta
use search, only: binary_search
@ -354,154 +352,136 @@ contains
integer, intent(in) :: i_nuclide ! index into nuclides array
real(8), intent(in) :: E ! energy
integer :: i ! loop index
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
integer :: same_nuc_idx ! index of same nuclide
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer :: urr
type(NuclideCE), pointer :: nuc
micro_xs(i_nuclide) % use_ptable = .true.
! get pointer to probability table
nuc => nuclides(i_nuclide)
urr => nuc % urr_data
associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data)
! determine energy table
i_energy = 1
do
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine energy table
i_energy = 1
do
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine interpolation factor on table
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! determine interpolation factor on table
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! sample probability table using the cumulative distribution
! sample probability table using the cumulative distribution
! Random numbers for xs calculation are sampled from a separated stream.
! This guarantees the randomness and, at the same time, makes sure we reuse
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
call prn_set_stream(STREAM_URR_PTABLE)
r = future_prn(int(nuc_zaid_dict % get_key(nuc % zaid), 8))
call prn_set_stream(STREAM_TRACKING)
! if we're dealing with a nuclide that we've previously encountered at
! this energy but a different temperature, use the original random number to
! preserve correlation of temperature in probability tables
same_nuc = .false.
do i = 1, nuc % nuc_list % size()
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % data(i)) % last_E) then
same_nuc = .true.
same_nuc_idx = i
exit
end if
end do
i_low = 1
do
if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
i_low = i_low + 1
end do
i_up = 1
do
if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
i_up = i_up + 1
end do
if (same_nuc) then
r = micro_xs(nuc % nuc_list % data(same_nuc_idx)) % last_prn
else
r = prn()
micro_xs(i_nuclide) % last_prn = r
end if
! determine elastic, fission, and capture cross sections from probability
! table
if (urr % interp == LINEAR_LINEAR) then
elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
f * urr % prob(i_energy + 1, URR_FISSION, i_up)
capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
i_low = 1
do
if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
i_low = i_low + 1
end do
i_up = 1
do
if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
i_up = i_up + 1
end do
! determine elastic, fission, and capture cross sections from probability
! table
if (urr % interp == LINEAR_LINEAR) then
elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
f * urr % prob(i_energy + 1, URR_FISSION, i_up)
capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
end if
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
end if
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
end if
end if
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine inelastic scattering cross section
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
end if
end associate
end if
! Multiply by smooth cross-section if needed
if (urr % multiply_smooth) then
elastic = elastic * micro_xs(i_nuclide) % elastic
capture = capture * (micro_xs(i_nuclide) % absorption - &
micro_xs(i_nuclide) % fission)
fission = fission * micro_xs(i_nuclide) % fission
end if
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
end if
! Check for negative values
if (elastic < ZERO) elastic = ZERO
if (fission < ZERO) fission = ZERO
if (capture < ZERO) capture = ZERO
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
end if
end if
! Set elastic, absorption, fission, and total cross sections. Note that the
! total cross section is calculated as sum of partials rather than using the
! table-provided value
micro_xs(i_nuclide) % elastic = elastic
micro_xs(i_nuclide) % absorption = capture + fission
micro_xs(i_nuclide) % fission = fission
micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine nu-fission cross section
if (nuc % fissionable) then
micro_xs(i_nuclide) % nu_fission = nu_total(nuc, E) * &
micro_xs(i_nuclide) % fission
end if
! Determine inelastic scattering cross section
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
end associate
end if
! Multiply by smooth cross-section if needed
if (urr % multiply_smooth) then
elastic = elastic * micro_xs(i_nuclide) % elastic
capture = capture * (micro_xs(i_nuclide) % absorption - &
micro_xs(i_nuclide) % fission)
fission = fission * micro_xs(i_nuclide) % fission
end if
! Check for negative values
if (elastic < ZERO) elastic = ZERO
if (fission < ZERO) fission = ZERO
if (capture < ZERO) capture = ZERO
! Set elastic, absorption, fission, and total cross sections. Note that the
! total cross section is calculated as sum of partials rather than using the
! table-provided value
micro_xs(i_nuclide) % elastic = elastic
micro_xs(i_nuclide) % absorption = capture + fission
micro_xs(i_nuclide) % fission = fission
micro_xs(i_nuclide) % total = elastic + inelastic + capture + fission
! Determine nu-fission cross section
if (nuc % fissionable) then
micro_xs(i_nuclide) % nu_fission = nuc % nu(E, EMISSION_TOTAL) * &
micro_xs(i_nuclide) % fission
end if
end associate
end subroutine calculate_urr_xs

4
src/eigenvalue.F90
View file

@ -10,7 +10,7 @@ module eigenvalue
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh
use random_lcg, only: prn, set_particle_seed, prn_skip
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use search, only: binary_search
use string, only: to_str
@ -99,7 +99,7 @@ contains
call set_particle_seed(int((current_batch - 1)*gen_per_batch + &
current_gen,8))
call prn_skip(start)
call advance_prn_seed(start)
! Determine how many fission sites we need to sample from the source bank
! and the probability for selecting a site.

190
src/endf_header.F90
View file

@ -1,12 +1,51 @@
module endf_header
implicit none
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG
use search, only: binary_search
implicit none
type, abstract :: Function1D
contains
procedure(function1d_evaluate_), deferred :: evaluate
end type Function1D
abstract interface
pure function function1d_evaluate_(this, x) result(y)
import Function1D
class(Function1D), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
end function function1d_evaluate_
end interface
!===============================================================================
! TAB1 represents a one-dimensional interpolable function
! CONSTANT1D represents a constant one-dimensional function
!===============================================================================
type Tab1
type, extends(Function1D) :: Constant1D
real(8) :: y
contains
procedure :: evaluate => constant1d_evaluate
end type Constant1D
!===============================================================================
! POLYNOMIAL represents a one-dimensional function expressed as a polynomial
!===============================================================================
type, extends(Function1D) :: Polynomial
real(8), allocatable :: coef(:) ! coefficients
contains
procedure :: evaluate => polynomial_evaluate
procedure :: from_ace => polynomial_from_ace
end type Polynomial
!===============================================================================
! TABULATED1D represents a one-dimensional interpolable function
!===============================================================================
type, extends(Function1D) :: Tabulated1D
integer :: n_regions = 0 ! # of interpolation regions
integer, allocatable :: nbt(:) ! values separating interpolation regions
integer, allocatable :: int(:) ! interpolation scheme
@ -14,18 +53,78 @@ module endf_header
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace
end type Tab1
procedure :: from_ace => tabulated1d_from_ace
procedure :: evaluate => tabulated1d_evaluate
end type Tabulated1D
contains
subroutine from_ace(this, xss, idx)
class(Tab1), intent(inout) :: this
!===============================================================================
! Constant1D implementation
!===============================================================================
pure function constant1d_evaluate(this, x) result(y)
class(Constant1D), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
y = this % y
end function constant1d_evaluate
!===============================================================================
! Polynomial implementation
!===============================================================================
subroutine polynomial_from_ace(this, xss, idx)
class(Polynomial), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nc ! number of coefficients (order - 1)
! Clear space
if (allocated(this % coef)) deallocate(this % coef)
! Determine number of coefficients
nc = nint(xss(idx))
! Allocate space for and read coefficients
allocate(this % coef(nc))
this % coef(:) = xss(idx + 1 : idx + nc)
end subroutine polynomial_from_ace
pure function polynomial_evaluate(this, x) result(y)
class(Polynomial), intent(in) :: this
real(8), intent(in) :: x
real(8) :: y
integer :: i
! Use Horner's rule to evaluate polynomial. Note that coefficients are
! ordered in increasing powers of x.
y = ZERO
do i = size(this % coef), 1, -1
y = y*x + this % coef(i)
end do
end function polynomial_evaluate
!===============================================================================
! Tabulated1D implementation
!===============================================================================
subroutine tabulated1d_from_ace(this, xss, idx)
class(Tabulated1D), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Clear space
if (allocated(this % nbt)) deallocate(this % nbt)
if (allocated(this % int)) deallocate(this % int)
if (allocated(this % x)) deallocate(this % x)
if (allocated(this % y)) deallocate(this % y)
! Determine number of regions
nr = nint(xss(idx))
this%n_regions = nr
@ -47,6 +146,81 @@ contains
allocate(this%y(ne))
this%x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne)
this%y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne)
end subroutine from_ace
end subroutine tabulated1d_from_ace
pure function tabulated1d_evaluate(this, x) result(y)
class(Tabulated1D), intent(in) :: this
real(8), intent(in) :: x ! x value to find y at
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine number of interpolation regions and pairs
n_regions = this % n_regions
n_pairs = this % n_pairs
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < this % x(1)) then
y = this % y(1)
return
elseif (x > this % x(n_pairs)) then
y = this % y(n_pairs)
return
else
i = binary_search(this % x, n_pairs, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = this % int(1)
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < this % nbt(j)) then
interp = this % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = this % y(i)
return
end if
! determine bounding values
x0 = this % x(i)
x1 = this % x(i + 1)
y0 = this % y(i)
y1 = this % y(i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function tabulated1d_evaluate
end module endf_header

125
src/energy_distribution.F90
View file

@ -1,8 +1,7 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use endf_header, only: Tabulated1D
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
@ -83,8 +82,8 @@ module energy_distribution
integer :: n_region
integer, allocatable :: breakpoints(:)
integer, allocatable :: interpolation(:)
real(8), allocatable :: energy_in(:)
type(CTTable), allocatable :: energy_out(:)
real(8), allocatable :: energy(:)
type(CTTable), allocatable :: distribution(:)
contains
procedure :: sample => continuous_sample
end type ContinuousTabular
@ -95,7 +94,7 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: MaxwellEnergy
type(Tab1) :: theta ! incoming-energy-dependent parameter
type(Tabulated1D) :: theta ! incoming-energy-dependent parameter
real(8) :: u ! restriction energy
contains
procedure :: sample => maxwellenergy_sample
@ -107,7 +106,7 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: Evaporation
type(Tab1) :: theta
type(Tabulated1D) :: theta
real(8) :: u
contains
procedure :: sample => evaporation_sample
@ -119,28 +118,13 @@ module energy_distribution
!===============================================================================
type, extends(EnergyDistribution) :: WattEnergy
type(Tab1) :: a
type(Tab1) :: b
type(Tabulated1D) :: a
type(Tabulated1D) :: b
real(8) :: u
contains
procedure :: sample => watt_sample
end type WattEnergy
!===============================================================================
! NBODYPHASESPACE gives the energy distribution for particles emitted from
! neutron and charged-particle reactions. This corresponds to ACE law 66 and
! ENDF File 6, LAW=6.
!===============================================================================
type, extends(EnergyDistribution) :: NBodyPhaseSpace
integer :: n_bodies
real(8) :: mass_ratio
real(8) :: A
real(8) :: Q
contains
procedure :: sample => nbody_sample
end type NBodyPhaseSpace
contains
function equiprobable_sample(this, E_in) result(E_out)
@ -202,6 +186,7 @@ contains
end if
end function equiprobable_sample
function level_inelastic_sample(this, E_in) result(E_out)
class(LevelInelastic), intent(in) :: this
real(8), intent(in) :: E_in
@ -210,6 +195,7 @@ contains
E_out = this%mass_ratio*(E_in - this%threshold)
end function level_inelastic_sample
function continuous_sample(this, E_in) result(E_out)
class(ContinuousTabular), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
@ -238,17 +224,17 @@ contains
! Find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
i = binary_search(this%energy, n_energy_in, E_in)
r = (E_in - this%energy(i)) / &
(this%energy(i+1) - this%energy(i))
end if
! Sample between the ith and (i+1)th bin
@ -263,23 +249,23 @@ contains
end if
! Interpolation for energy E1 and EK
n_energy_out = size(this%energy_out(i)%e_out)
E_i_1 = this%energy_out(i)%e_out(1)
E_i_K = this%energy_out(i)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i)%e_out)
E_i_1 = this%distribution(i)%e_out(1)
E_i_K = this%distribution(i)%e_out(n_energy_out)
n_energy_out = size(this%energy_out(i+1)%e_out)
E_i1_1 = this%energy_out(i+1)%e_out(1)
E_i1_K = this%energy_out(i+1)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i+1)%e_out)
E_i1_1 = this%distribution(i+1)%e_out(1)
E_i1_K = this%distribution(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! Determine outgoing energy bin
n_energy_out = size(this%energy_out(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%energy_out(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%energy_out(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -287,9 +273,9 @@ contains
! Check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%energy_out(l)%e_out(k)
p_l_k = this%energy_out(l)%p(k)
if (this%energy_out(l)%interpolation == HISTOGRAM) then
E_l_k = this%distribution(l)%e_out(k)
p_l_k = this%distribution(l)%p(k)
if (this%distribution(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
@ -297,10 +283,10 @@ contains
E_out = E_l_k
end if
elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%energy_out(l)%e_out(k+1)
p_l_k1 = this%energy_out(l)%p(k+1)
E_l_k1 = this%distribution(l)%e_out(k+1)
p_l_k1 = this%distribution(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
@ -321,6 +307,7 @@ contains
end if
end function continuous_sample
function maxwellenergy_sample(this, E_in) result(E_out)
class(MaxwellEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
@ -329,7 +316,7 @@ contains
real(8) :: theta ! Maxwell distribution parameter
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
theta = this % theta % evaluate(E_in)
do
! Sample maxwell fission spectrum
@ -349,9 +336,9 @@ contains
real(8) :: x, y, v
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
theta = this % theta % evaluate(E_in)
y = (E_in - this%U)/theta
y = (E_in - this%u)/theta
v = 1 - exp(-y)
! Sample outgoing energy based on evaporation spectrum probability
@ -372,10 +359,10 @@ contains
real(8) :: a, b ! Watt spectrum parameters
! Determine Watt parameter 'a' from tabulated function
a = interpolate_tab1(this%a, E_in)
a = this % a % evaluate(E_in)
! Determine Watt parameter 'b' from tabulated function
b = interpolate_tab1(this%b, E_in)
b = this % b % evaluate(E_in)
do
! Sample energy-dependent Watt fission spectrum
@ -386,44 +373,4 @@ contains
end do
end function watt_sample
function nbody_sample(this, E_in) result(E_out)
class(NBodyPhaseSpace), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: Ap ! total mass of particles in neutron masses
real(8) :: E_max ! maximum possible COM energy
real(8) :: x, y, v
real(8) :: r1, r2, r3, r4, r5, r6
! Determine E_max parameter
Ap = this%mass_ratio
E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q)
! x is essentially a Maxwellian distribution
x = maxwell_spectrum(ONE)
select case (this%n_bodies)
case (3)
y = maxwell_spectrum(ONE)
case (4)
r1 = prn()
r2 = prn()
r3 = prn()
y = -log(r1*r2*r3)
case (5)
r1 = prn()
r2 = prn()
r3 = prn()
r4 = prn()
r5 = prn()
r6 = prn()
y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
end select
! Now determine v and E_out
v = x/(x+y)
E_out = E_max * v
end function nbody_sample
end module energy_distribution

161
src/fission.F90
View file

@ -1,161 +0,0 @@
module fission
use nuclide_header, only: NuclideCE
use constants
use error, only: fatal_error
use interpolation, only: interpolate_tab1
use search, only: binary_search
implicit none
contains
!===============================================================================
! NU_TOTAL calculates the total number of neutrons emitted per fission for a
! given nuclide and incoming neutron energy
!===============================================================================
pure function nu_total(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_t_type == NU_NONE) then
nu = ERROR_REAL
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_t_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_t_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_t_data, E)
end if
end function nu_total
!===============================================================================
! NU_PROMPT calculates the total number of prompt neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
pure function nu_prompt(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_p_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call nu_total unnecessarily if it has already been called.
nu = ZERO
elseif (nuc % nu_p_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_p_data(1))
! sum up polynomial in energy
nu = ZERO
do i = 0, NC - 1
c = nuc % nu_p_data(i+2)
nu = nu + c * E**i
end do
elseif (nuc % nu_p_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_p_data, E)
end if
end function nu_prompt
!===============================================================================
! NU_DELAYED calculates the total number of delayed neutrons emitted per fission
! for a given nuclide and incoming neutron energy
!===============================================================================
pure function nu_delayed(nuc, E) result(nu)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
if (nuc % nu_d_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call nu_delayed unnecessarily if it has already been called.
nu = ZERO
elseif (nuc % nu_d_type == NU_TABULAR) then
! use ENDF interpolation laws to determine nu
nu = interpolate_tab1(nuc % nu_d_data, E)
end if
end function nu_delayed
!===============================================================================
! YIELD_DELAYED calculates the fractional yield of delayed neutrons emitted for
! a given nuclide and incoming neutron energy in a given delayed group.
!===============================================================================
pure function yield_delayed(nuc, E, g) result(yield)
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: yield ! delayed neutron precursor yield
integer, intent(in) :: g ! the delayed neutron precursor group
integer :: d ! precursor group
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
yield = ZERO
if (g > nuc % n_precursor .or. g < 1) then
! if the precursor group is outside the range of precursor groups for
! the input nuclide, return ZERO.
yield = ZERO
else if (nuc % nu_d_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
! emission is prompt -- WARNING: This currently returns zero. The calling
! routine needs to know this situation is occurring since we don't want
! to call yield_delayed unnecessarily if it has already been called.
yield = ZERO
else if (nuc % nu_d_type == NU_TABULAR) then
lc = 1
! loop over delayed groups and determine the yield for the desired group
do d = 1, nuc % n_precursor
! determine number of interpolation regions and energies
NR = int(nuc % nu_d_precursor_data(lc + 1))
NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR))
! check if this is the desired group
if (d == g) then
! determine delayed neutron precursor yield for group g
yield = interpolate_tab1(nuc % nu_d_precursor_data( &
lc+1:lc+2+2*NR+2*NE), E)
exit
end if
! advance pointer
lc = lc + 2 + 2*NR + 2*NE + 1
end do
end if
end function yield_delayed
end module fission

4
src/global.F90
View file

@ -105,6 +105,10 @@ module global
! What to assume for expanding natural elements
integer :: default_expand = ENDF_BVII1
! Total amount of nuclide ZAID and dictionary of nuclide ZAID and index
integer(8) :: n_nuc_zaid_total
type(DictIntInt) :: nuc_zaid_dict
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES

11
src/initialize.F90
View file

@ -1,6 +1,6 @@
module initialize
use ace, only: read_ace_xs, same_nuclide_list
use ace, only: read_ace_xs
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
@ -16,7 +16,7 @@ module initialize
hdf5_tallyresult_t, hdf5_integer8_t
use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use mgxs_data, only: read_mgxs, same_NuclideMG_list, create_macro_xs
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: title, header, print_version, write_message, &
print_usage, write_xs_summary, print_plot
use random_lcg, only: initialize_prng
@ -122,13 +122,6 @@ contains
end if
call time_read_xs%stop()
! Create linked lists for multiple instances of the same nuclide
if (run_CE) then
call same_nuclide_list()
else
call same_nuclidemg_list()
end if
! Construct information needed for nuclear data
if (run_CE) then
! Construct unionized or log energy grid for cross-sections

105
src/input_xml.F90
View file

@ -1891,21 +1891,23 @@ contains
subroutine read_materials_xml()
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index for elements
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: n_nuc_ele ! number of nuclides in an element
integer :: index_list ! index in xs_listings array
integer :: index_nuclide ! index in nuclides
integer :: index_sab ! index in sab_tables
real(8) :: val ! value entered for density
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index for elements
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: n_nuc_ele ! number of nuclides in an element
integer :: index_list ! index in xs_listings array
integer :: index_nuclide ! index in nuclides
integer :: index_nuc_zaid ! index in nuclide ZAID
integer :: index_sab ! index in sab_tables
real(8) :: val ! value entered for density
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
integer :: zaid ! ZAID of nuclide
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
character(MAX_WORD_LEN) :: units ! units on density
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
@ -1955,6 +1957,7 @@ contains
! Initialize count for number of nuclides/S(a,b) tables
index_nuclide = 0
index_nuc_zaid = 0
index_sab = 0
do i = 1, n_materials
@ -2095,21 +2098,6 @@ contains
end if
end if
! Check enforced isotropic lab scattering
if (check_for_node(node_nuc, "scattering")) then
call get_node_value(node_nuc, "scattering", temp_str)
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
call list_iso_lab % append(1)
else if (adjustl(to_lower(temp_str)) == "data") then
call list_iso_lab % append(0)
else
call fatal_error("Scattering must be isotropic in lab or follow&
& the ACE file data")
end if
else
call list_iso_lab % append(0)
end if
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
@ -2154,6 +2142,23 @@ contains
end if
end if
! Check enforced isotropic lab scattering
if (run_CE) then
if (check_for_node(node_nuc, "scattering")) then
call get_node_value(node_nuc, "scattering", temp_str)
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
call list_iso_lab % append(1)
else if (adjustl(to_lower(temp_str)) == "data") then
call list_iso_lab % append(0)
else
call fatal_error("Scattering must be isotropic in lab or follow&
& the ACE file data")
end if
else
call list_iso_lab % append(0)
end if
end if
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
@ -2248,23 +2253,25 @@ contains
n_nuc_ele = list_names % size() - n_nuc_ele
! Check enforced isotropic lab scattering
if (check_for_node(node_ele, "scattering")) then
call get_node_value(node_ele, "scattering", temp_str)
else
temp_str = "data"
end if
! Set ace or iso-in-lab scattering for each nuclide in element
do k = 1, n_nuc_ele
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
call list_iso_lab % append(1)
else if (adjustl(to_lower(temp_str)) == "data") then
call list_iso_lab % append(0)
if (run_CE) then
if (check_for_node(node_ele, "scattering")) then
call get_node_value(node_ele, "scattering", temp_str)
else
call fatal_error("Scattering must be isotropic in lab or follow&
& the ACE file data")
temp_str = "data"
end if
end do
! Set ace or iso-in-lab scattering for each nuclide in element
do k = 1, n_nuc_ele
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
call list_iso_lab % append(1)
else if (adjustl(to_lower(temp_str)) == "data") then
call list_iso_lab % append(0)
else
call fatal_error("Scattering must be isotropic in lab or follow&
& the ACE file data")
end if
end do
end if
end do NATURAL_ELEMENTS
@ -2300,6 +2307,7 @@ contains
index_list = xs_listing_dict % get_key(to_lower(name))
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
zaid = xs_listings(index_list) % zaid
! If this nuclide hasn't been encountered yet, we need to add its name
! and alias to the nuclide_dict
@ -2313,6 +2321,12 @@ contains
mat % nuclide(j) = nuclide_dict % get_key(to_lower(name))
end if
! Construct dict of nuclide zaid
if (.not. nuc_zaid_dict % has_key(zaid)) then
index_nuc_zaid = index_nuc_zaid + 1
call nuc_zaid_dict % add_key(zaid, index_nuc_zaid)
end if
! Copy name and atom/weight percent
mat % names(j) = name
mat % atom_density(j) = list_density % get_item(j)
@ -2407,6 +2421,7 @@ contains
! Set total number of nuclides and S(a,b) tables
n_nuclides_total = index_nuclide
n_sab_tables = index_sab
n_nuc_zaid_total = index_nuc_zaid
! Close materials XML file
call close_xmldoc(doc)

205
src/interpolation.F90
View file

@ -1,205 +0,0 @@
module interpolation
use constants
use endf_header, only: Tab1
use search, only: binary_search
use string, only: to_str
implicit none
interface interpolate_tab1
module procedure interpolate_tab1_array, interpolate_tab1_object
end interface interpolate_tab1
contains
!===============================================================================
! INTERPOLATE_TAB1_ARRAY interpolates a function between two points based on
! particular interpolation scheme. The data needs to be organized as a ENDF TAB1
! type function containing the interpolation regions, break points, and
! tabulated x's and y's.
!===============================================================================
pure function interpolate_tab1_array(data, x, loc_start) result(y)
real(8), intent(in) :: data(:) ! array of data
real(8), intent(in) :: x ! x value to find y at
integer, intent(in), optional :: loc_start ! starting location in data
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: loc_0 ! starting location
integer :: n_regions ! number of interpolation regions
integer :: n_points ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
integer :: loc_breakpoints ! location of breakpoints in data
integer :: loc_interp ! location of interpolation schemes in data
integer :: loc_x ! location of x's in data
integer :: loc_y ! location of y's in data
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine starting location
if (present(loc_start)) then
loc_0 = loc_start - 1
else
loc_0 = 0
end if
! determine number of interpolation regions
n_regions = int(data(loc_0 + 1))
! set locations for breakpoints and interpolation schemes
loc_breakpoints = loc_0 + 1
loc_interp = loc_breakpoints + n_regions
! determine number of tabulated values
n_points = int(data(loc_interp + n_regions + 1))
! set locations for x's and y's
loc_x = loc_interp + n_regions + 1
loc_y = loc_x + n_points
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < data(loc_x + 1)) then
y = data(loc_y + 1)
return
elseif (x > data(loc_x + n_points)) then
y = data(loc_y + n_points)
return
else
i = binary_search(data(loc_x + 1:loc_x + n_points), n_points, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = int(data(loc_interp + 1))
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < data(loc_breakpoints + j)) then
interp = int(data(loc_interp + j))
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = data(loc_y + i)
return
end if
! determine bounding values
x0 = data(loc_x + i)
x1 = data(loc_x + i + 1)
y0 = data(loc_y + i)
y1 = data(loc_y + i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_array
!===============================================================================
! INTERPOLATE_TAB1_OBJECT interpolates a function between two points based on
! particular interpolation scheme. The data needs to be organized as a ENDF TAB1
! type function containing the interpolation regions, break points, and
! tabulated x's and y's.
!===============================================================================
pure function interpolate_tab1_object(obj, x) result(y)
type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function
real(8), intent(in) :: x ! x value to find y at
real(8) :: y ! y(x)
integer :: i ! bin in which to interpolate
integer :: j ! index for interpolation region
integer :: n_regions ! number of interpolation regions
integer :: n_pairs ! number of tabulated values
integer :: interp ! ENDF interpolation scheme
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
! determine number of interpolation regions and pairs
n_regions = obj % n_regions
n_pairs = obj % n_pairs
! find which bin the abscissa is in -- if the abscissa is outside the
! tabulated range, the first or last point is chosen, i.e. no interpolation
! is done outside the energy range
if (x < obj % x(1)) then
y = obj % y(1)
return
elseif (x > obj % x(n_pairs)) then
y = obj % y(n_pairs)
return
else
i = binary_search(obj % x, n_pairs, x)
end if
! determine interpolation scheme
if (n_regions == 0) then
interp = LINEAR_LINEAR
elseif (n_regions == 1) then
interp = obj % int(1)
elseif (n_regions > 1) then
do j = 1, n_regions
if (i < obj % nbt(j)) then
interp = obj % int(j)
exit
end if
end do
end if
! handle special case of histogram interpolation
if (interp == HISTOGRAM) then
y = obj % y(i)
return
end if
! determine bounding values
x0 = obj % x(i)
x1 = obj % x(i + 1)
y0 = obj % y(i)
y1 = obj % y(i + 1)
! determine interpolation factor and interpolated value
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_object
end module interpolation

22
src/mgxs_data.F90
View file

@ -161,28 +161,6 @@ contains
end subroutine read_mgxs
!===============================================================================
! SAME_NUCLIDEMG_LIST creates a linked list for each nuclide containing the
! indices in the nuclides array of all other instances of that nuclide. For
! example, the same nuclide may exist at multiple temperatures resulting
! in multiple entries in the nuclides array for a single zaid number.
!===============================================================================
subroutine same_nuclidemg_list()
integer :: i ! index in nuclides array
integer :: j ! index in nuclides array
do i = 1, n_nuclides_total
do j = 1, n_nuclides_total
if (nuclides_MG(i) % obj % zaid == nuclides_MG(j) % obj % zaid) then
call nuclides_MG(i) % obj % nuc_list % push_back(j)
end if
end do
end do
end subroutine same_nuclidemg_list
!===============================================================================
! CREATE_MACRO_XS generates the macroscopic x/s from the microscopic input data
!===============================================================================

216
src/nuclide_header.F90
View file

@ -2,13 +2,18 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header
use constants
use endf, only: reaction_name
use error, only: fatal_error
use dict_header, only: DictIntInt
use endf, only: reaction_name, is_fission, is_disappearance
use endf_header, only: Function1D
use error, only: fatal_error, warning
use list_header, only: ListInt
use math, only: evaluate_legendre, find_angle
use product_header, only: AngleEnergyContainer
use reaction_header, only: Reaction
use stl_vector, only: VectorInt
use string
use urr_header, only: UrrData
use xml_interface
implicit none
@ -26,9 +31,6 @@ module nuclide_header
integer :: listing ! index in xs_listings
real(8) :: kT ! temperature in MeV (k*T)
! Linked list of indices in nuclides array of instances of this same nuclide
type(VectorInt) :: nuc_list
! Fission information
logical :: fissionable ! nuclide is fissionable?
@ -70,24 +72,11 @@ module nuclide_header
real(8) :: E_max ! upper cutoff energy for res scattering
! Fission information
logical :: has_partial_fission ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
logical :: has_partial_fission = .false. ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
integer :: n_precursor = 0 ! # of delayed neutron precursors
integer, allocatable :: index_fission(:) ! indices in reactions
! Total fission neutron emission
integer :: nu_t_type
real(8), allocatable :: nu_t_data(:)
! Prompt fission neutron emission
integer :: nu_p_type
real(8), allocatable :: nu_p_data(:)
! Delayed fission neutron emission
integer :: nu_d_type
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
class(Function1D), allocatable :: total_nu
! Unresolved resonance data
logical :: urr_present
@ -103,6 +92,7 @@ module nuclide_header
contains
procedure :: clear => nuclidece_clear
procedure :: print => nuclidece_print
procedure :: nu => nuclidece_nu
end type NuclideCE
type, abstract, extends(Nuclide) :: NuclideMG
@ -257,7 +247,6 @@ module nuclide_header
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
real(8) :: last_prn
end type NuclideMicroXS
!===============================================================================
@ -684,81 +673,133 @@ module nuclide_header
! or NuclideAngle
!===============================================================================
subroutine nuclidece_clear(this)
subroutine nuclidece_clear(this)
class(NuclideCE), intent(inout) :: this ! The Nuclide object to clear
class(NuclideCE), intent(inout) :: this ! The Nuclide object to clear
integer :: i ! Loop counter
if (associated(this % urr_data)) deallocate(this % urr_data)
if (associated(this % urr_data)) deallocate(this % urr_data)
call this % reaction_index % clear()
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()
end do
end subroutine nuclidece_clear
function nuclidece_nu(this, E, emission_mode, group) result(nu)
class(NuclideCE), intent(in) :: this
real(8), intent(in) :: E
integer, intent(in) :: emission_mode
integer, optional, intent(in) :: group
real(8) :: nu
integer :: i
if (.not. this % fissionable) then
nu = ZERO
return
end if
select case (emission_mode)
case (EMISSION_PROMPT)
associate (product => this % reactions(this % index_fission(1)) % products(1))
nu = product % yield % evaluate(E)
end associate
case (EMISSION_DELAYED)
if (this % n_precursor > 0) then
if (present(group)) then
! If delayed group specified, determine yield immediately
associate(p => this % reactions(this % index_fission(1)) % products(1 + group))
nu = p % yield % evaluate(E)
end associate
else
nu = ZERO
associate (rx => this % reactions(this % index_fission(1)))
do i = 2, size(rx % products)
associate (product => rx % products(i))
! Skip any non-neutron products
if (product % particle /= NEUTRON) exit
! Evaluate yield
if (product % emission_mode == EMISSION_DELAYED) then
nu = nu + product % yield % evaluate(E)
end if
end associate
end do
end associate
end if
else
nu = ZERO
end if
call this % reaction_index % clear()
case (EMISSION_TOTAL)
if (allocated(this % total_nu)) then
nu = this % total_nu % evaluate(E)
else
associate (rx => this % reactions(this % index_fission(1)))
nu = rx % products(1) % yield % evaluate(E)
end associate
end if
end select
end subroutine nuclidece_clear
end function nuclidece_nu
!===============================================================================
! NUCLIDE*_PRINT displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine nuclidece_print(this, unit)
class(NuclideCE), intent(in) :: this
integer, intent(in), optional :: unit
subroutine nuclidece_print(this, unit)
class(NuclideCE), intent(in) :: this
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
type(UrrData), pointer :: urr
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Initialize totals
size_urr = 0
size_xs = 0
! Initialize totals
size_urr = 0
size_xs = 0
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction))
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, this % n_reaction
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, this % n_reaction
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * this % n_grid * 4
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * this % n_grid * 4
! Write information about URR probability tables
size_urr = 0
if (this % urr_present) then
urr => this % urr_data
! Write information about URR probability tables
size_urr = 0
if (this % urr_present) then
associate(urr => this % urr_data)
write(unit_,*) ' Unresolved resonance probability table:'
write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy))
write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob))
@ -771,17 +812,18 @@ module nuclide_header
! Calculate memory used by probability tables and add to total
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end if
end associate
end if
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine nuclidece_print
! Blank line at end of nuclide
write(unit_,*)
end subroutine nuclidece_print
subroutine nuclidemg_print(this, unit_)
class(NuclideMG), intent(in) :: this

3
src/output.F90
View file

@ -2,7 +2,6 @@ module output
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning
@ -1176,7 +1175,7 @@ contains
function get_label(t, i_filter) result(label)
type(TallyObject), intent(in) :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(100) :: label ! user-specified identifier
character(MAX_LINE_LEN) :: label ! user-specified identifier
integer :: i ! index in cells/surfaces/etc array
integer :: bin

141
src/physics.F90
View file

@ -1,13 +1,10 @@
module physics
use ace_header, only: Reaction
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
use error, only: fatal_error, warning
use fission, only: nu_total, nu_delayed
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use math
use mesh, only: get_mesh_indices
@ -16,7 +13,8 @@ module physics
use particle_header, only: Particle
use particle_restart_write, only: write_particle_restart
use physics_common
use random_lcg, only: prn
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
use reaction_header, only: Reaction
use search, only: binary_search
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str
@ -58,6 +56,13 @@ contains
if (master) call warning("Killing neutron with extremely low energy")
end if
! Advance URR seed stream 'N' times after energy changes
if (p % E /= p % last_E) then
call prn_set_stream(STREAM_URR_PTABLE)
call advance_prn_seed(n_nuc_zaid_total)
call prn_set_stream(STREAM_TRACKING)
endif
end subroutine collision
!===============================================================================
@ -440,9 +445,9 @@ contains
vel = sqrt(dot_product(v_n, v_n))
! Sample scattering angle
select type (dist => rxn%secondary%distribution(1)%obj)
select type (dist => rxn % products(1) % distribution(1) % obj)
type is (UncorrelatedAngleEnergy)
mu_cm = dist%angle%sample(E)
mu_cm = dist % angle % sample(E)
end select
! Determine direction cosines in CM
@ -1066,8 +1071,6 @@ contains
integer :: nu ! actual number of neutrons produced
integer :: ijk(3) ! indices in ufs mesh
real(8) :: nu_t ! total nu
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(NuclideCE), pointer :: nuc
@ -1138,25 +1141,12 @@ contains
! Set weight of fission bank site
bank_array(i) % wgt = ONE/weight
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Sample delayed group and angle/energy for fission reaction
call sample_fission_neutron(nuc, nuc % reactions(i_reaction), &
p % E, bank_array(i))
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
bank_array(i) % uvw(1) = mu
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
bank_array(i) % E = sample_fission_energy(nuc, &
nuc % reactions(i_reaction), p)
! Set the delayed group of the neutron
bank_array(i) % delayed_group = p % delayed_group
! Set delayed group on particle too
p % delayed_group = bank_array(i) % delayed_group
! Increment the number of neutrons born delayed
if (p % delayed_group > 0) then
@ -1175,35 +1165,41 @@ contains
end subroutine create_fission_sites
!===============================================================================
! SAMPLE_FISSION_ENERGY
! SAMPLE_FISSION_NEUTRON
!===============================================================================
function sample_fission_energy(nuc, rxn, p) result(E_out)
subroutine sample_fission_neutron(nuc, rxn, E_in, site)
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
real(8), intent(in) :: E_in
type(Bank), intent(inout) :: site
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing energy of fission neutron
integer :: j ! index on nu energy grid / precursor group
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
integer :: n_sample ! number of times resampling
integer :: group ! index on nu energy grid / precursor group
integer :: n_sample ! number of resamples
real(8) :: nu_t ! total nu
real(8) :: nu_d ! delayed nu
real(8) :: beta ! delayed neutron fraction
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
real(8) :: mu ! cosine of scattering angle
real(8) :: phi ! azimuthal angle
! Determine total nu
nu_t = nu_total(nuc, p % E)
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Determine delayed nu
nu_d = nu_delayed(nuc, p % E)
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
site % uvw(1) = mu
site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Determine delayed neutron fraction
! Determine total nu, delayed nu, and delayed neutron fraction
nu_t = nuc % nu(E_in, EMISSION_TOTAL)
nu_d = nuc % nu(E_in, EMISSION_DELAYED)
beta = nu_d / nu_t
if (prn() < beta) then
@ -1211,51 +1207,41 @@ contains
! DELAYED NEUTRON SAMPLED
! sampled delayed precursor group
xi = prn()
lc = 1
xi = prn()*nu_d
prob = ZERO
do j = 1, nuc % n_precursor
! determine number of interpolation regions and energies
NR = int(nuc % nu_d_precursor_data(lc + 1))
NE = int(nuc % nu_d_precursor_data(lc + 2 + 2*NR))
do group = 1, nuc % n_precursor
! determine delayed neutron precursor yield for group j
yield = interpolate_tab1(nuc % nu_d_precursor_data( &
lc+1:lc+2+2*NR+2*NE), p % E)
yield = rxn % products(1 + group) % yield % evaluate(E_in)
! Check if this group is sampled
prob = prob + yield
if (xi < prob) exit
! advance pointer
lc = lc + 2 + 2*NR + 2*NE + 1
end do
! if the sum of the probabilities is slightly less than one and the
! random number is greater, j will be greater than nuc %
! n_precursor -- check for this condition
j = min(j, nuc % n_precursor)
group = min(group, nuc % n_precursor)
! set the delayed group for the particle born from fission
p % delayed_group = j
site % delayed_group = group
! sample from energy distribution
n_sample = 0
do
select type (aedist => nuc%nu_d_edist(j)%obj)
type is (UncorrelatedAngleEnergy)
E_out = aedist%energy%sample(p%E)
end select
! sample from energy/angle distribution -- note that mu has already been
! sampled above and doesn't need to be resampled
call rxn % products(1 + group) % sample(E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
if (site % E < energy_max_neutron) exit
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call write_particle_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
&// "times for nuclide " // nuc % name)
// "times for nuclide " // nuc % name)
end if
end do
@ -1264,28 +1250,27 @@ contains
! PROMPT NEUTRON SAMPLED
! set the delayed group for the particle born from fission to 0
p % delayed_group = 0
site % delayed_group = 0
! sample from prompt neutron energy distribution
n_sample = 0
do
call rxn%secondary%sample(p%E, E_out, prob)
call rxn % products(1) % sample(E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
if (site % E < energy_max_neutron) exit
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call write_particle_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
&// "times for nuclide " // nuc % name)
// "times for nuclide " // nuc % name)
end if
end do
end if
end function sample_fission_energy
end subroutine sample_fission_neutron
!===============================================================================
! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other
@ -1309,7 +1294,7 @@ contains
E_in = p % E
! sample outgoing energy and scattering cosine
call rxn%secondary%sample(E_in, E, mu)
call rxn % products(1) % sample(E_in, E, mu)
! if scattering system is in center-of-mass, transfer cosine of scattering
! angle and outgoing energy from CM to LAB
@ -1337,14 +1322,16 @@ contains
! change direction of particle
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
! change weight of particle based on yield
if (rxn % multiplicity_with_E) then
yield = interpolate_tab1(rxn % multiplicity_E, E_in)
p % wgt = yield * p % wgt
else
do i = 1, rxn % multiplicity - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.True.)
! evaluate yield
yield = rxn % products(1) % yield % evaluate(E_in)
if (mod(yield, ONE) == ZERO) then
! If yield is integral, create exactly that many secondary particles
do i = 1, nint(yield) - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.true.)
end do
else
! Otherwise, change weight of particle based on yield
p % wgt = yield * p % wgt
end if
end subroutine inelastic_scatter

62
src/product_header.F90 Normal file
View file

@ -0,0 +1,62 @@
module product_header
use angleenergy_header, only: AngleEnergyContainer
use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, &
EMISSION_TOTAL, NEUTRON, PHOTON
use endf_header, only: Tabulated1D, Function1D, Constant1D, Polynomial
use random_lcg, only: prn
!===============================================================================
! REACTIONPRODUCT stores a data for a reaction product including its yield and
! angle-energy distributions, each of which has a given probability of occurring
! for a given incoming energy. In general, most products only have one
! angle-energy distribution, but for some cases (e.g., (n,2n) in certain
! nuclides) multiple distinct distributions exist.
!===============================================================================
type :: ReactionProduct
integer :: particle
integer :: emission_mode ! prompt, delayed, or total emission
real(8) :: decay_rate ! Decay rate for delayed neutron precursors
class(Function1D), pointer :: yield => null() ! Energy-dependent neutron yield
type(Tabulated1D), allocatable :: applicability(:)
type(AngleEnergyContainer), allocatable :: distribution(:)
contains
procedure :: sample => reactionproduct_sample
end type ReactionProduct
contains
subroutine reactionproduct_sample(this, E_in, E_out, mu)
class(ReactionProduct), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
integer :: i ! loop counter
integer :: n ! number of angle-energy distributions
real(8) :: prob ! cumulative probability
real(8) :: c ! sampled cumulative probability
n = size(this%applicability)
if (n > 1) then
prob = ZERO
c = prn()
do i = 1, n
! Determine probability that i-th energy distribution is sampled
prob = prob + this % applicability(i) % evaluate(E_in)
! If i-th distribution is sampled, sample energy from the distribution
if (c <= prob) then
call this%distribution(i)%obj%sample(E_in, E_out, mu)
exit
end if
end do
else
! If only one distribution is present, go ahead and sample it
call this%distribution(1)%obj%sample(E_in, E_out, mu)
end if
end subroutine reactionproduct_sample
end module product_header

35
src/random_lcg.F90
View file

@ -24,9 +24,10 @@ module random_lcg
!$omp threadprivate(prn_seed, stream)
public :: prn
public :: future_prn
public :: initialize_prng
public :: set_particle_seed
public :: prn_skip
public :: advance_prn_seed
public :: prn_set_stream
public :: STREAM_TRACKING, STREAM_TALLIES
@ -52,6 +53,21 @@ contains
end function prn
!===============================================================================
! FUTURE_PRN generates a pseudo-random number which is 'n' times ahead from the
! current seed.
!===============================================================================
function future_prn(n) result(pseudo_rn)
integer(8), intent(in) :: n ! number of prns to skip
real(8) :: pseudo_rn
pseudo_rn = future_seed(n, prn_seed(stream)) * prn_norm
end function future_prn
!===============================================================================
! INITIALIZE_PRNG sets up the random number generator, determining the seed and
! values for g, c, and m.
@ -90,31 +106,32 @@ contains
integer :: i
do i = 1, N_STREAMS
prn_seed(i) = prn_skip_ahead(id*prn_stride, prn_seed0 + i - 1)
prn_seed(i) = future_seed(id*prn_stride, prn_seed0 + i - 1)
end do
end subroutine set_particle_seed
!===============================================================================
! PRN_SKIP advances the random number seed 'n' times from the current seed
! ADVANCE_PRN_SEED advances the random number seed 'n' times from the current
! seed.
!===============================================================================
subroutine prn_skip(n)
subroutine advance_prn_seed(n)
integer(8), intent(in) :: n ! number of seeds to skip
prn_seed(stream) = prn_skip_ahead(n, prn_seed(stream))
prn_seed(stream) = future_seed(n, prn_seed(stream))
end subroutine prn_skip
end subroutine advance_prn_seed
!===============================================================================
! PRN_SKIP_AHEAD advances the random number seed 'skip' times. This is usually
! FUTURE_SEED advances the random number seed 'skip' times. This is usually
! used to skip a fixed number of random numbers (the stride) so that a given
! particle always has the same starting seed regardless of how many processors
! are used
!===============================================================================
function prn_skip_ahead(n, seed) result(new_seed)
function future_seed(n, seed) result(new_seed)
integer(8), intent(in) :: n ! number of seeds to skip
integer(8), intent(in) :: seed ! original seed
@ -166,7 +183,7 @@ contains
! With G and C, we can now find the new seed
new_seed = iand(g_new*seed + c_new, prn_mask)
end function prn_skip_ahead
end function future_seed
!===============================================================================
! PRN_SET_STREAM changes the random number stream. If random numbers are needed

21
src/reaction_header.F90 Normal file
View file

@ -0,0 +1,21 @@
module reaction_header
use product_header, only: ReactionProduct
implicit none
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
real(8), allocatable :: sigma(:) ! Cross section values
type(ReactionProduct), allocatable :: products(:)
end type Reaction
end module reaction_header

52
src/secondary_correlated.F90
View file

@ -1,8 +1,8 @@
module secondary_correlated
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
@ -24,8 +24,8 @@ module secondary_correlated
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(AngleEnergyTable), allocatable :: table(:) ! outgoing E/mu distributions
real(8), allocatable :: energy(:) ! incoming energies
type(AngleEnergyTable), allocatable :: distribution(:) ! outgoing E/mu distributions
contains
procedure :: sample => correlated_sample
end type CorrelatedAngleEnergy
@ -61,17 +61,17 @@ contains
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
i = binary_search(this%energy, n_energy_in, E_in)
r = (E_in - this%energy(i)) / &
(this%energy(i+1) - this%energy(i))
end if
! Sample between the ith and (i+1)th bin
@ -82,23 +82,23 @@ contains
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i)%e_out)
E_i_1 = this%distribution(i)%e_out(1)
E_i_K = this%distribution(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i+1)%e_out)
E_i1_1 = this%distribution(i+1)%e_out(1)
E_i1_K = this%distribution(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -106,9 +106,9 @@ contains
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
E_l_k = this%distribution(l)%e_out(k)
p_l_k = this%distribution(l)%p(k)
if (this%distribution(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
@ -116,10 +116,10 @@ contains
E_out = E_l_k
end if
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
E_l_k1 = this%distribution(l)%e_out(k+1)
p_l_k1 = this%distribution(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
@ -139,9 +139,9 @@ contains
! Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1) then
mu = this%table(l)%angle(k)%obj%sample()
mu = this%distribution(l)%angle(k)%obj%sample()
else
mu = this%table(l)%angle(k + 1)%obj%sample()
mu = this%distribution(l)%angle(k + 1)%obj%sample()
end if
end subroutine correlated_sample

83
src/secondary_header.F90
View file

@ -1,83 +0,0 @@
module secondary_header
use constants, only: ZERO
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use random_lcg, only: prn
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
! distribution that is a function of incoming energy. Each derived type must
! implement a sample() subroutine that returns an outgoing energy and scattering
! cosine given an incoming energy.
!===============================================================================
type, abstract :: AngleEnergy
contains
procedure(angleenergy_sample_), deferred :: sample
end type AngleEnergy
abstract interface
subroutine angleenergy_sample_(this, E_in, E_out, mu)
import AngleEnergy
class(AngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in
real(8), intent(out) :: E_out
real(8), intent(out) :: mu
end subroutine angleenergy_sample_
end interface
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
!===============================================================================
! SECONDARYDISTRIBUTION stores multiple angle-energy distributions, each of
! which has a given probability of occurring for a given incoming energy. In
! general, most secondary distributions only have one angle-energy distribution,
! but for some cases (e.g., (n,2n) in certain nuclides) multiple distinct
! distributions exist.
!===============================================================================
type :: SecondaryDistribution
type(Tab1), allocatable :: applicability(:)
type(AngleEnergyContainer), allocatable :: distribution(:)
contains
procedure :: sample => secondary_sample
end type SecondaryDistribution
contains
subroutine secondary_sample(this, E_in, E_out, mu)
class(SecondaryDistribution), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
integer :: i ! loop counter
integer :: n ! number of angle-energy distributions
real(8) :: prob ! cumulative probability
real(8) :: c ! sampled cumulative probability
n = size(this%applicability)
if (n > 1) then
prob = ZERO
c = prn()
do i = 1, n
! Determine probability that i-th energy distribution is sampled
prob = prob + interpolate_tab1(this%applicability(i), E_in)
! If i-th distribution is sampled, sample energy from the distribution
if (c <= prob) then
call this%distribution(i)%obj%sample(E_in, E_out, mu)
exit
end if
end do
else
! If only one distribution is present, go ahead and sample it
call this%distribution(1)%obj%sample(E_in, E_out, mu)
end if
end subroutine secondary_sample
end module secondary_header

60
src/secondary_kalbach.F90
View file

@ -1,7 +1,7 @@
module secondary_kalbach
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
@ -25,8 +25,8 @@ module secondary_kalbach
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(KalbachMannTable), allocatable :: table(:) ! outgoing E/mu parameters
real(8), allocatable :: energy(:) ! incoming energies
type(KalbachMannTable), allocatable :: distribution(:) ! outgoing E/mu parameters
contains
procedure :: sample => kalbachmann_sample
end type KalbachMann
@ -64,17 +64,17 @@ contains
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
n_energy_in = size(this%energy)
if (E_in < this%energy(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
elseif (E_in > this%energy(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
i = binary_search(this%energy, n_energy_in, E_in)
r = (E_in - this%energy(i)) / &
(this%energy(i+1) - this%energy(i))
end if
! Sample between the ith and (i+1)th bin
@ -85,23 +85,23 @@ contains
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i)%e_out)
E_i_1 = this%distribution(i)%e_out(1)
E_i_K = this%distribution(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
n_energy_out = size(this%distribution(i+1)%e_out)
E_i1_1 = this%distribution(i+1)%e_out(1)
E_i1_K = this%distribution(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
n_energy_out = size(this%distribution(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
c_k = this%distribution(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
c_k1 = this%distribution(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
@ -109,9 +109,9 @@ contains
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
E_l_k = this%distribution(l)%e_out(k)
p_l_k = this%distribution(l)%p(k)
if (this%distribution(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
@ -120,13 +120,13 @@ contains
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k)
km_a = this%table(l)%a(k)
km_r = this%distribution(l)%r(k)
km_a = this%distribution(l)%a(k)
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
E_l_k1 = this%distribution(l)%e_out(k+1)
p_l_k1 = this%distribution(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
@ -137,10 +137,10 @@ contains
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%r(k+1) - this%table(l)%r(k))
km_a = this%table(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%a(k+1) - this%table(l)%a(k))
km_r = this%distribution(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%distribution(l)%r(k+1) - this%distribution(l)%r(k))
km_a = this%distribution(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%distribution(l)%a(k+1) - this%distribution(l)%a(k))
end if
! Now interpolate between incident energy bins i and i + 1

70
src/secondary_nbody.F90 Normal file
View file

@ -0,0 +1,70 @@
module secondary_nbody
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, PI
use math, only: maxwell_spectrum
use random_lcg, only: prn
!===============================================================================
! NBODYPHASESPACE gives the energy distribution for particles emitted from
! neutron and charged-particle reactions. This corresponds to ACE law 66 and
! ENDF File 6, LAW=6.
!===============================================================================
type, extends(AngleEnergy) :: NBodyPhaseSpace
integer :: n_bodies
real(8) :: mass_ratio
real(8) :: A
real(8) :: Q
contains
procedure :: sample => nbody_sample
end type NBodyPhaseSpace
contains
subroutine nbody_sample(this, E_in, E_out, mu)
class(NBodyPhaseSpace), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled outgoing energy
real(8) :: Ap ! total mass of particles in neutron masses
real(8) :: E_max ! maximum possible COM energy
real(8) :: x, y, v
real(8) :: r1, r2, r3, r4, r5, r6
! By definition, the distribution of the angle is isotropic for an N-body
! phase space distribution
mu = TWO*prn() - ONE
! Determine E_max parameter
Ap = this%mass_ratio
E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q)
! x is essentially a Maxwellian distribution
x = maxwell_spectrum(ONE)
select case (this%n_bodies)
case (3)
y = maxwell_spectrum(ONE)
case (4)
r1 = prn()
r2 = prn()
r3 = prn()
y = -log(r1*r2*r3)
case (5)
r1 = prn()
r2 = prn()
r3 = prn()
r4 = prn()
r5 = prn()
r6 = prn()
y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
end select
! Now determine v and E_out
v = x/(x+y)
E_out = E_max * v
end subroutine nbody_sample
end module secondary_nbody

2
src/secondary_uncorrelated.F90
View file

@ -1,9 +1,9 @@
module secondary_uncorrelated
use angle_distribution, only: AngleDistribution
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO
use energy_distribution, only: EnergyDistribution
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
!===============================================================================

148
src/state_point.F90
View file

@ -49,10 +49,9 @@ contains
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group, tally_group
integer(HID_T) :: meshes_group, mesh_group
integer(HID_T) :: filter_group, derivs_group, deriv_group
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
mesh_group, filter_group, derivs_group, deriv_group, &
runtime_group
character(20), allocatable :: str_array(:)
character(MAX_FILE_LEN) :: filename
type(RegularMesh), pointer :: meshp
@ -133,13 +132,13 @@ contains
call write_dataset(file_id, "cmfd_on", 1)
cmfd_group = create_group(file_id, "cmfd")
call write_dataset(cmfd_group, "indices", cmfd%indices)
call write_dataset(cmfd_group, "k_cmfd", cmfd%k_cmfd)
call write_dataset(cmfd_group, "cmfd_src", cmfd%cmfd_src)
call write_dataset(cmfd_group, "cmfd_entropy", cmfd%entropy)
call write_dataset(cmfd_group, "cmfd_balance", cmfd%balance)
call write_dataset(cmfd_group, "cmfd_dominance", cmfd%dom)
call write_dataset(cmfd_group, "cmfd_srccmp", cmfd%src_cmp)
call write_dataset(cmfd_group, "indices", cmfd % indices)
call write_dataset(cmfd_group, "k_cmfd", cmfd % k_cmfd)
call write_dataset(cmfd_group, "cmfd_src", cmfd % cmfd_src)
call write_dataset(cmfd_group, "cmfd_entropy", cmfd % entropy)
call write_dataset(cmfd_group, "cmfd_balance", cmfd % balance)
call write_dataset(cmfd_group, "cmfd_dominance", cmfd % dom)
call write_dataset(cmfd_group, "cmfd_srccmp", cmfd % src_cmp)
call close_group(cmfd_group)
else
call write_dataset(file_id, "cmfd_on", 0)
@ -155,18 +154,18 @@ contains
if (n_meshes > 0) then
! Print list of mesh IDs
current => mesh_dict%keys()
current => mesh_dict % keys()
allocate(id_array(n_meshes))
allocate(key_array(n_meshes))
i = 1
do while (associated(current))
key_array(i) = current%key
id_array(i) = current%value
key_array(i) = current % key
id_array(i) = current % value
! Move to next mesh
next => current%next
next => current % next
deallocate(current)
current => next
i = i + 1
@ -180,16 +179,17 @@ contains
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
meshp => meshes(id_array(i))
mesh_group = create_group(meshes_group, "mesh " // trim(to_str(meshp%id)))
mesh_group = create_group(meshes_group, "mesh " &
// trim(to_str(meshp % id)))
select case (meshp%type)
select case (meshp % type)
case (MESH_REGULAR)
call write_dataset(mesh_group, "type", "regular")
end select
call write_dataset(mesh_group, "dimension", meshp%dimension)
call write_dataset(mesh_group, "lower_left", meshp%lower_left)
call write_dataset(mesh_group, "upper_right", meshp%upper_right)
call write_dataset(mesh_group, "width", meshp%width)
call write_dataset(mesh_group, "dimension", meshp % dimension)
call write_dataset(mesh_group, "lower_left", meshp % lower_left)
call write_dataset(mesh_group, "upper_right", meshp % upper_right)
call write_dataset(mesh_group, "width", meshp % width)
call close_group(mesh_group)
end do MESH_LOOP
@ -240,7 +240,7 @@ contains
! Write all tally information except results
do i = 1, n_tallies
tally => tallies(i)
key_array(i) = tally%id
key_array(i) = tally % id
id_array(i) = i
end do
@ -255,9 +255,9 @@ contains
! Get pointer to tally
tally => tallies(i)
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally%id)))
trim(to_str(tally % id)))
select case(tally%estimator)
select case(tally % estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
case (ESTIMATOR_TRACKLENGTH)
@ -265,16 +265,17 @@ contains
case (ESTIMATOR_COLLISION)
call write_dataset(tally_group, "estimator", "collision")
end select
call write_dataset(tally_group, "n_realizations", tally%n_realizations)
call write_dataset(tally_group, "n_filters", tally%n_filters)
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
call write_dataset(tally_group, "n_filters", tally % n_filters)
! Write filter information
FILTER_LOOP: do j = 1, tally%n_filters
FILTER_LOOP: do j = 1, tally % n_filters
filter_group = create_group(tally_group, "filter " // &
trim(to_str(j)))
! Write name of type
select case (tally%filters(j)%type)
select case (tally % filters(j) % type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
@ -303,36 +304,37 @@ contains
call write_dataset(filter_group, "type", "delayedgroup")
end select
call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins)
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. &
tally % filters(j) % type == FILTER_MU .or. &
tally % filters(j) % type == FILTER_POLAR .or. &
tally % filters(j) % type == FILTER_AZIMUTHAL) then
call write_dataset(filter_group, "bins", &
tally%filters(j)%real_bins)
tally % filters(j) % real_bins)
else
call write_dataset(filter_group, "bins", &
tally%filters(j)%int_bins)
tally % filters(j) % int_bins)
end if
call close_group(filter_group)
end do FILTER_LOOP
! Set up nuclide bin array and then write
allocate(str_array(tally%n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally%n_nuclide_bins
if (tally%nuclide_bins(j) > 0) then
allocate(str_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (tally % nuclide_bins(j) > 0) then
! Get index in cross section listings for this nuclide
i_list = nuclides(tally%nuclide_bins(j))%listing
i_list = nuclides(tally % nuclide_bins(j)) % listing
! Determine position of . in alias string (e.g. "U-235.71c"). If
! no . is found, just use the entire string.
i_xs = index(xs_listings(i_list)%alias, '.')
i_xs = index(xs_listings(i_list) % alias, '.')
if (i_xs > 0) then
str_array(j) = xs_listings(i_list)%alias(1:i_xs - 1)
str_array(j) = xs_listings(i_list) % alias(1:i_xs - 1)
else
str_array(j) = xs_listings(i_list)%alias
str_array(j) = xs_listings(i_list) % alias
end if
else
str_array(j) = 'total'
@ -348,32 +350,33 @@ contains
end if
! Write scores.
call write_dataset(tally_group, "n_score_bins", tally%n_score_bins)
allocate(str_array(size(tally%score_bins)))
do j = 1, size(tally%score_bins)
str_array(j) = reaction_name(tally%score_bins(j))
call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
allocate(str_array(size(tally % score_bins)))
do j = 1, size(tally % score_bins)
str_array(j) = reaction_name(tally % score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
call write_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins)
call write_dataset(tally_group, "n_user_score_bins", &
tally % n_user_score_bins)
deallocate(str_array)
! Write explicit moment order strings for each score bin
k = 1
allocate(str_array(tally%n_score_bins))
MOMENT_LOOP: do j = 1, tally%n_user_score_bins
select case(tally%score_bins(k))
allocate(str_array(tally % n_score_bins))
MOMENT_LOOP: do j = 1, tally % n_user_score_bins
select case(tally % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
str_array(k) = 'P' // trim(to_str(tally%moment_order(k)))
str_array(k) = 'P' // trim(to_str(tally % moment_order(k)))
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, tally%moment_order(k)
do n_order = 0, tally % 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, tally%moment_order(k)
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
str_array(k) = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
@ -425,8 +428,9 @@ contains
tally => tallies(i)
! Write sum and sum_sq for each bin
tally_group = open_group(tallies_group, "tally " // to_str(tally%id))
call write_dataset(tally_group, "results", tally%results)
tally_group = open_group(tallies_group, "tally " &
// to_str(tally % id))
call write_dataset(tally_group, "results", tally % results)
call close_group(tally_group)
end do TALLY_RESULTS
@ -436,13 +440,45 @@ contains
end if
call close_group(tallies_group)
! Write out the runtime metrics.
runtime_group = create_group(file_id, "runtime")
call write_dataset(runtime_group, "total initialization", &
time_initialize % get_value())
call write_dataset(runtime_group, "reading cross sections", &
time_read_xs % get_value())
call write_dataset(runtime_group, "simulation", &
time_inactive % get_value() + time_active % get_value())
call write_dataset(runtime_group, "transport", &
time_transport % get_value())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "inactive batches", &
time_inactive % get_value())
end if
call write_dataset(runtime_group, "active batches", &
time_active % get_value())
if (run_mode == MODE_EIGENVALUE) then
call write_dataset(runtime_group, "synchronizing fission bank", &
time_bank % get_value())
call write_dataset(runtime_group, "sampling source sites", &
time_bank_sample % get_value())
call write_dataset(runtime_group, "SEND-RECV source sites", &
time_bank_sendrecv % get_value())
end if
call write_dataset(runtime_group, "accumulating tallies", &
time_tallies % get_value())
if (cmfd_run) then
call write_dataset(runtime_group, "CMFD", time_cmfd % get_value())
call write_dataset(runtime_group, "CMFD building matrices", &
time_cmfdbuild % get_value())
call write_dataset(runtime_group, "CMFD solving matrices", &
time_cmfdsolve % get_value())
end if
call write_dataset(runtime_group, "total", time_total % get_value())
call close_group(runtime_group)
call file_close(file_id)
end if
if (master .and. n_tallies > 0) then
deallocate(id_array)
end if
end subroutine write_state_point
!===============================================================================

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