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

Merge branch 'develop' into integrator-class

Browse source 
Need new travis.yml to pass CI
This commit is contained in:
Andrew Johnson 2019-07-30 15:58:49 -05:00
commit b9683dbdba
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
114 changed files with 8323 additions and 4462 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.travis.yml
View file

@ -28,6 +28,7 @@ env:
- LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib
- PATH=$PATH:$HOME/NJOY2016/build
- COVERALLS_PARALLEL=true
- NUMPY_EXPERIMENTAL_ARRAY_FUNCTION=0
matrix:
include:
- python: "3.5"

56
docs/source/devguide/tests.rst
View file

@ -4,9 +4,6 @@
Test Suite
==========
Running Tests
-------------
The OpenMC test suite consists of two parts, a regression test suite and a unit
test suite. The regression test suite is based on regression or integrated
testing where different types of input files are configured and the full OpenMC
@ -14,27 +11,35 @@ code is executed. Results from simulations are compared with expected
results. The unit tests are primarily intended to test individual
functions/classes in the OpenMC Python API.
The test suite relies on the third-party `pytest <https://pytest.org>`_
package. To run either or both the regression and unit test suites, it is
assumed that you have OpenMC fully installed, i.e., the :ref:`scripts_openmc`
executable is available on your :envvar:`PATH` and the :mod:`openmc` Python
module is importable. In development where it would be onerous to continually
install OpenMC every time a small change is made, it is recommended to install
OpenMC in development/editable mode. With setuptools, this is accomplished by
running::
Prerequisites
-------------
python setup.py develop
- The test suite relies on the third-party `pytest <https://pytest.org>`_
package. To run either or both the regression and unit test suites, it is
assumed that you have OpenMC fully installed, i.e., the :ref:`scripts_openmc`
executable is available on your :envvar:`PATH` and the :mod:`openmc` Python
module is importable. In development where it would be onerous to continually
install OpenMC every time a small change is made, it is recommended to install
OpenMC in development/editable mode. With setuptools, this is accomplished by
running::
or using pip (recommended)::
python setup.py develop
pip install -e .[test]
or using pip (recommended)::
It is also assumed that you have cross section data available that is pointed to
by the :envvar:`OPENMC_CROSS_SECTIONS` environment variables. Furthermore, to
run unit tests for the :mod:`openmc.data` module, it is necessary to have
ENDF/B-VII.1 data available and pointed to by the :envvar:`OPENMC_ENDF_DATA`
environment variable. All data sources can be obtained using the
``tools/ci/travis-before-script.sh`` script.
pip install -e .[test]
- The test suite requires a specific set of cross section data in order for
tests to pass. A download URL for the data that OpenMC expects can be found
within ``tools/ci/download-xs.sh``.
- In addition to the HDF5 data, some tests rely on ENDF files. A download URL
for those can also be found in ``tools/ci/download-xs.sh``.
- Some tests require `NJOY <https://www.njoy21.io/NJOY2016>`_ to preprocess
cross section data. The test suite assumes that you have an ``njoy``
executable available on your :envvar:`PATH`.
Running Tests
-------------
To execute the test suite, go to the ``tests/`` directory and run::
@ -46,6 +51,17 @@ installed and run::
pytest --cov=../openmc --cov-report=html
Generating XML Inputs
---------------------
Many of the regression tests rely on the Python API to build an appropriate
model. However, it can sometimes be desirable to work directly with the XML
input files rather than having to run a script in order to run the problem/test.
To build the input files for a test without actually running the test, you can
run::
pytest --build-inputs <name-of-test>
Adding Tests to the Regression Suite
------------------------------------

17
docs/source/io_formats/cross_sections.rst
View file

@ -51,3 +51,20 @@ attributes:
:type:
The type of data contained in the file. Accepted values are 'neutron',
'thermal', 'photon', and 'wmp'.
.. _depletion_element:
-----------------------------
``<depletion_chain>`` Element
-----------------------------
The ``<depletion_chain>`` element indicates the location of the depletion chain file.
This file contains information describing how nuclides decay and transmute to other
nuclides through the depletion process. This element has a single attribute, ``path``,
pointing to the location of the chain file.
.. code-block:: xml
<depletion_chain path="/opt/data/chain_endfb7.xml"/>
The structure of the depletion chain file is explained in :ref:`io_depletion_chain`.

7
docs/source/io_formats/depletion_results.rst
View file

@ -44,3 +44,10 @@ The current version of the depletion results file format is 1.0.
**/reactions/<name>/**
:Attributes: - **index** (*int*) -- Index user in results for this reaction
.. note::
The reaction rates for some isotopes not originally present may
be non-zero, but should be negligible compared to other atoms.
This can be controlled by changing the
:class:`openmc.deplete.Operator` ``dilute_initial`` attribute.

3
docs/source/io_formats/materials.rst
View file

@ -42,8 +42,7 @@ Each ``material`` element can have the following attributes or sub-elements:
Volume of the material in cm^3.
:temperature:
An element with no attributes which is used to set the default temperature
of the material in Kelvin.
Temperature of the material in Kelvin.
*Default*: If a material default temperature is not given and a cell
temperature is not specified, the :ref:`global default temperature

7
docs/source/io_formats/source.rst
View file

@ -13,8 +13,9 @@ is that documented here.
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
:Datasets:
- **source_bank** (Compound type) -- Source bank information for each
particle. The compound type has fields ``wgt``, ``xyz``, ``uvw``,
``E``, and ``delayed_group``, which represent the weight, position,
direction, energy, energy group, and delayed_group of the source
particle, respectively.
``E``, ``delayed_group``, and ``particle``, which represent the
weight, position, direction, energy, energy group, delayed group,
and type of the source particle, respectively.

7
docs/source/io_formats/summary.rst
View file

@ -116,6 +116,13 @@ The current version of the summary file format is 6.0.
- **sab_names** (*char[][]*) -- Names of
S(:math:`\alpha,\beta`) tables assigned to the material.
:Attributes: - **volume** (*double[]*) -- Volume of this material [cm^3]. Only
present if ``volume`` supplied
- **temperature** (*double[]*) -- Temperature of this material [K].
Only present in ``temperature`` supplied
- **depletable** (*int[]*) -- ``1`` if the material can be depleted,
``0`` otherwise. Always present
**/nuclides/**
:Attributes: - **n_nuclides** (*int*) -- Number of nuclides in the problem.

14
docs/source/pythonapi/deplete.rst
View file

@ -75,10 +75,24 @@ data, such as number densities and reaction rates for each material.
:template: myclass.rst
AtomNumber
ChainFissionHelper
DirectReactionRateHelper
OperatorResult
ReactionRates
Results
ResultsList
The following classes are abstract classes that can be used to extend the
:mod:`openmc.deplete` capabilities:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
ReactionRateHelper
EnergyHelper
TransportOperator
Each of the integrator functions also relies on a number of "helper" functions

1
docs/source/pythonapi/model.rst
View file

@ -15,6 +15,7 @@ Convenience Functions
openmc.model.hexagonal_prism
openmc.model.rectangular_prism
openmc.model.subdivide
openmc.model.pin
TRISO Fuel Modeling
-------------------

7
examples/jupyter/cad-based-geometry.ipynb
View file

@ -236,6 +236,13 @@
"tallies.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Note:** Applying tally filters in DagMC models requires prior knowledge of the model. Here, we know that the fuel cell's volume ID in the CAD sofware is 1. To identify cells without use of CAD software, load them into the [OpenMC plotter](https://github.com/openmc/plotter) where cell, material, and volume IDs can be identified for native both OpenMC and DagMC geometries."
]
},
{
"cell_type": "markdown",
"metadata": {},

697
examples/jupyter/candu.ipynb
View file

File diff suppressed because one or more lines are too long

550
examples/jupyter/mdgxs-part-i.ipynb
View file

File diff suppressed because one or more lines are too long

460
examples/jupyter/mdgxs-part-ii.ipynb
View file

File diff suppressed because one or more lines are too long

80
examples/jupyter/mgxs-part-i.ipynb
View file

@ -28,9 +28,7 @@
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -134,9 +132,7 @@
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
@ -224,9 +220,7 @@
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Instantiate a Cell\n",
@ -268,9 +262,7 @@
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Create Geometry and set root Universe\n",
@ -326,9 +318,7 @@
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Instantiate a 2-group EnergyGroups object\n",
@ -365,9 +355,7 @@
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Instantiate a few different sections\n",
@ -390,9 +378,7 @@
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -431,9 +417,7 @@
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"name": "stderr",
@ -473,9 +457,7 @@
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"name": "stdout",
@ -645,9 +627,7 @@
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Load the last statepoint file\n",
@ -671,9 +651,7 @@
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"# Load the tallies from the statepoint into each MGXS object\n",
@ -706,9 +684,7 @@
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"name": "stdout",
@ -741,9 +717,7 @@
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -807,9 +781,7 @@
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"absorption.export_xs_data(filename='absorption-xs', format='excel')"
@ -825,9 +797,7 @@
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [],
"source": [
"total.build_hdf5_store(filename='mgxs', append=True)\n",
@ -852,9 +822,7 @@
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -931,9 +899,7 @@
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -1003,9 +969,7 @@
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -1082,9 +1046,7 @@
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"metadata": {},
"outputs": [
{
"data": {
@ -1168,9 +1130,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.0"
"version": "3.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
"nbformat_minor": 1
}

2412
examples/jupyter/mgxs-part-ii.ipynb
View file

File diff suppressed because one or more lines are too long

871
examples/jupyter/mgxs-part-iii.ipynb
View file

File diff suppressed because one or more lines are too long

162
examples/jupyter/nuclear-data-resonance-covariance.ipynb
View file

File diff suppressed because one or more lines are too long

142
examples/jupyter/nuclear-data.ipynb
View file

File diff suppressed because one or more lines are too long

1233
examples/jupyter/pandas-dataframes.ipynb
View file

File diff suppressed because one or more lines are too long

557
examples/jupyter/pincell.ipynb
View file

File diff suppressed because one or more lines are too long

515
examples/jupyter/post-processing.ipynb
View file

File diff suppressed because one or more lines are too long

26
examples/jupyter/search.ipynb
View file

@ -35,7 +35,7 @@
"\n",
"To perform the search we will use the `openmc.search_for_keff` function. This function requires a different function be defined which creates an parametrized model to analyze. This model is required to be stored in an `openmc.model.Model` object. The first parameter of this function will be modified during the search process for our critical eigenvalue.\n",
"\n",
"Our model will be a pin-cell from the [Multi-Group Mode Part II](http://openmc.readthedocs.io/en/latest/examples/mg-mode-part-ii.html) assembly, except this time the entire model building process will be contained within a function, and the Boron concentration will be parametrized."
"Our model will be a pin-cell from the [Multi-Group Mode Part II](http://docs.openmc.org/en/latest/examples/mg-mode-part-ii.html) assembly, except this time the entire model building process will be contained within a function, and the Boron concentration will be parametrized."
]
},
{
@ -64,7 +64,7 @@
"\n",
" # Include the amount of boron in the water based on the ppm,\n",
" # neglecting the other constituents of boric acid\n",
" water.add_element('B', ppm_Boron * 1E-6)\n",
" water.add_element('B', ppm_Boron * 1e-6)\n",
" \n",
" # Instantiate a Materials object\n",
" materials = openmc.Materials([fuel, zircaloy, water])\n",
@ -141,7 +141,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -156,7 +156,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -171,7 +171,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -186,7 +186,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -201,7 +201,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -216,7 +216,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -231,7 +231,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -246,7 +246,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -261,7 +261,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/home/shriwise/.pyenv/versions/3.6.7/lib/python3.6/site-packages/openmc-0.10.0-py3.6-linux-x86_64.egg/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Universe instance already exists with id=0.\n",
" warn(msg, IDWarning)\n"
]
},
@ -277,7 +277,7 @@
"source": [
"# Perform the search\n",
"crit_ppm, guesses, keffs = openmc.search_for_keff(build_model, bracket=[1000., 2500.],\n",
" tol=1.E-2, bracketed_method='bisect',\n",
" tol=1e-2, bracketed_method='bisect',\n",
" print_iterations=True)\n",
"\n",
"print('Critical Boron Concentration: {:4.0f} ppm'.format(crit_ppm))"
@ -344,7 +344,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
"version": "3.7.0"
}
},
"nbformat": 4,

498
examples/jupyter/tally-arithmetic.ipynb
View file

@ -10,9 +10,7 @@
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"import glob\n",
@ -39,9 +37,7 @@
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# 1.6 enriched fuel\n",
@ -74,9 +70,7 @@
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Instantiate a Materials collection\n",
@ -96,14 +90,12 @@
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Create cylinders for the fuel and clad\n",
"fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218)\n",
"clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720)\n",
"fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.39218)\n",
"clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.45720)\n",
"\n",
"# Create boundary planes to surround the geometry\n",
"# Use both reflective and vacuum boundaries to make life interesting\n",
@ -125,9 +117,7 @@
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Create a Universe to encapsulate a fuel pin\n",
@ -162,9 +152,7 @@
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Create root Cell\n",
@ -189,9 +177,7 @@
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Create Geometry and set root Universe\n",
@ -201,9 +187,7 @@
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Export to \"geometry.xml\"\n",
@ -220,9 +204,7 @@
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# OpenMC simulation parameters\n",
@ -256,10 +238,19 @@
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEVyEhLpgJFNv8T///98iBL0AAAAAWJLR0QDEQxM8gAAAAd0SU1FB+MHEhEzAhO1TdMAAAKlSURBVGje7ZrBscIwDETxwSWkn5TAgXCgBPqhBA6kyj/fDhCIJa2zZAwz0pmHpZWdSazd7Tw8PDw8PDw8vinCMBzW03HIsRLvhnv0HL7qD+IwjzXKzaNaxeEt9kz21RWEBV5XQbfka3pQWL4qgdLyNQkUcbwFT/FP4zjWt+D+++OY4laZQJgtnqNOwe5l9XkGmIIL/PEHUAGTeuc5P15wBbu34ucSIAXkX77h4xUtIBSXnxIAOhCLy98TANNfLj8lYBYQCuLPWmAWEBe9f90DUPdKy08JWB0U1HsoaAgQxPSnAgwBopz+VABQvoDnAnqTP0r8zealzfPcQqqAQSs/C6AKGLX0pwKs8uX0cwGaAHr6uYC9wSt46qDCB4RXBDTkMwWMhnxZQF3+s8pf1AZY5VsCWuVnAUQ+YPxB43X5koAiH035soCa/AaeBOw34m359AaQPCK/1oAAyJ8aIPBI+7QGRkD+3IBt+A6QPzeg34SH2pcauN+Kt9uXGljkse0jb6BP8AD+vwGKPLZ95A0UofbnDbAFj20/eQN+gD8h/LgRD25/8QCA2088AD/Oo8dPOoDo8ZMOoPPNeej4pwdAgUcfX9IDzHnnf5lnz88XnH/nSf4M8cIL7I+/P3yCP0G88P7W+v2z9ft36+8P9vuJ/X5r/f3Jfj83//5vff/R+v6Hvb9i78/Y+7vW94/N71/Z+2P2/pq9P2fv7+n5ATu/YOcn7PyGnR+x8yt6ftYN3PzOENCcH7LzS3Z+Ss9vO62DV5uPmgAXSz5+fs7O72n/QBQLwPwLrH+C9W/Q/hHWv8L6Z2j/ThZgvX+I9S/R/inWv8X6x2j/Guufo/17rH+Q9S/S/knWv0n7R2n/Kuufpf27tH+Y9i/vWP+0h4eHh4eHh8cW8QcxLJDBvLKoigAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOS0wNy0xOFQyMjo1MTowMi0wNTowMAMmdtQAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTktMDctMThUMjI6NTE6MDItMDU6MDBye85oAAAAAElFTkSuQmCC\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Instantiate a Plot\n",
"plot = openmc.Plot(plot_id=1)\n",
@ -269,51 +260,8 @@
"plot.pixels = [250, 250]\n",
"plot.color_by = 'material'\n",
"\n",
"# Instantiate a Plots collection and export to \"plots.xml\"\n",
"plot_file = openmc.Plots([plot])\n",
"plot_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the plots.xml file, we can now generate and view the plot. OpenMC outputs plots in .ppm format, which can be converted into a compressed format like .png with the convert utility."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Run openmc in plotting mode\n",
"openmc.plot_geometry(output=False)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAxQTFRF\n////chIS6YCRTb/E6kGE+wAAAAFiS0dEAIgFHUgAAAAJcEhZcwAAAEgAAABIAEbJaz4AAALKSURB\nVGje7dpLcqQwDAbgHHE2YeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmN\nP+HDhw8fPnz48Kf6VH9G+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4\nzPji99z0/AJ4n1lfvJ6fnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6\npA0wfln+ho/fwgYYn19C/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tN\nDbSGz7T0SBEWw4vLXzbQ6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X5\n8wZaxWd1+fMGiuFvir8bvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV\n873hB8UnM3xzANtf8nb4dwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7\nT/ppARBvp48UwJnelT5SACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4/\n/Jve+fhsH6Ctv7n8PTzjvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V\n32/o9+fl389Xnx+g5x/o+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6\n/4Le/6D3T/D9V67Y/ZsVQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/\ngPs/0P4TtP8F7r9J3AIO9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTu\nf4X7b+H+X7T/+BPuf3aM8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTgtMDQtMDNUMjE6MTE6MzgtMDQ6MDD1dVTHAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE4LTA0LTAz\nVDIxOjExOjM4LTA0OjAwhCjsewAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Convert OpenMC's funky ppm to png\n",
"!convert materials-xy.ppm materials-xy.png\n",
"\n",
"# Display the materials plot inline\n",
"Image(filename='materials-xy.png')"
"# Show plot\n",
"openmc.plot_inline(plot)"
]
},
{
@ -325,10 +273,8 @@
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Instantiate an empty Tallies object\n",
@ -337,10 +283,8 @@
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Create Tallies to compute microscopic multi-group cross-sections\n",
@ -397,10 +341,8 @@
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# K-Eigenvalue (infinity) tallies\n",
@ -413,10 +355,8 @@
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": true
},
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Resonance Escape Probability tallies\n",
@ -428,10 +368,8 @@
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# Thermal Flux Utilization tallies\n",
@ -444,10 +382,8 @@
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": true
},
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# Fast Fission Factor tallies\n",
@ -459,10 +395,8 @@
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true
},
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# Instantiate energy filter to illustrate Tally slicing\n",
@ -479,18 +413,18 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/liangjg/.local/lib/python3.5/site-packages/openmc-0.10.0-py3.5.egg/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=6.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=6.\n",
" warn(msg, IDWarning)\n",
"/home/liangjg/.local/lib/python3.5/site-packages/openmc-0.10.0-py3.5.egg/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=3.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=3.\n",
" warn(msg, IDWarning)\n",
"/home/liangjg/.local/lib/python3.5/site-packages/openmc-0.10.0-py3.5.egg/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=2.\n",
"/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Filter instance already exists with id=2.\n",
" warn(msg, IDWarning)\n"
]
}
@ -509,7 +443,7 @@
},
{
"cell_type": "code",
"execution_count": 21,
"execution_count": 19,
"metadata": {
"scrolled": true
},
@ -518,66 +452,63 @@
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" %%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ################## %%%%%%%%%%%%%%%%%%%%%%%\n",
" ################### %%%%%%%%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%%%%%%\n",
" ##################### %%%%%%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%\n",
" ################# %%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%\n",
" ############ %%%%%%%%%%%%%%%\n",
" ######## %%%%%%%%%%%%%%\n",
" %%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ################## %%%%%%%%%%%%%%%%%%%%%%%\n",
" ################### %%%%%%%%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%%%%%%\n",
" ##################### %%%%%%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%\n",
" ################# %%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%\n",
" ############ %%%%%%%%%%%%%%%\n",
" ######## %%%%%%%%%%%%%%\n",
" %%%%%%%%%%%\n",
"\n",
" | The OpenMC Monte Carlo Code\n",
" Copyright | 2011-2018 Massachusetts Institute of Technology\n",
" Copyright | 2011-2019 MIT and OpenMC contributors\n",
" License | http://openmc.readthedocs.io/en/latest/license.html\n",
" Version | 0.10.0\n",
" Git SHA1 | 47fbf8282ea94c138f75219bd10fdb31501d3fb7\n",
" Date/Time | 2018-04-03 21:12:27\n",
" MPI Processes | 1\n",
" OpenMP Threads | 20\n",
" Version | 0.11.0-dev\n",
" Git SHA1 | 61c911cffdae2406f9f4bc667a9a6954748bb70c\n",
" Date/Time | 2019-07-18 22:51:02\n",
" OpenMP Threads | 4\n",
"\n",
" Reading settings XML file...\n",
" Reading cross sections XML file...\n",
" Reading materials XML file...\n",
" Reading geometry XML file...\n",
" Building neighboring cells lists for each surface...\n",
" Reading U235 from /home/liangjg/nucdata/nndc_hdf5/U235.h5\n",
" Reading U238 from /home/liangjg/nucdata/nndc_hdf5/U238.h5\n",
" Reading O16 from /home/liangjg/nucdata/nndc_hdf5/O16.h5\n",
" Reading H1 from /home/liangjg/nucdata/nndc_hdf5/H1.h5\n",
" Reading B10 from /home/liangjg/nucdata/nndc_hdf5/B10.h5\n",
" Reading Zr90 from /home/liangjg/nucdata/nndc_hdf5/Zr90.h5\n",
" Maximum neutron transport energy: 2.00000E+07 eV for U235\n",
" Reading U235 from /opt/data/hdf5/nndc_hdf5_v15/U235.h5\n",
" Reading U238 from /opt/data/hdf5/nndc_hdf5_v15/U238.h5\n",
" Reading O16 from /opt/data/hdf5/nndc_hdf5_v15/O16.h5\n",
" Reading H1 from /opt/data/hdf5/nndc_hdf5_v15/H1.h5\n",
" Reading B10 from /opt/data/hdf5/nndc_hdf5_v15/B10.h5\n",
" Reading Zr90 from /opt/data/hdf5/nndc_hdf5_v15/Zr90.h5\n",
" Maximum neutron transport energy: 20000000.000000 eV for U235\n",
" Reading tallies XML file...\n",
" Writing summary.h5 file...\n",
" Initializing source particles...\n",
"\n",
" ====================> K EIGENVALUE SIMULATION <====================\n",
"\n",
" Bat./Gen. k Average k \n",
" ========= ======== ==================== \n",
" 1/1 0.96168 \n",
" 2/1 0.96651 \n",
" 3/1 1.00678 \n",
" 4/1 0.98773 \n",
" 5/1 1.01883 \n",
" 6/1 1.02959 \n",
" Bat./Gen. k Average k\n",
" ========= ======== ====================\n",
" 1/1 0.96168\n",
" 2/1 0.96651\n",
" 3/1 1.00678\n",
" 4/1 0.98773\n",
" 5/1 1.01883\n",
" 6/1 1.02959\n",
" 7/1 0.99859 1.01409 +/- 0.01550\n",
" 8/1 1.03441 1.02086 +/- 0.01123\n",
" 9/1 1.06097 1.03089 +/- 0.01279\n",
@ -596,28 +527,28 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 4.9090E-01 seconds\n",
" Reading cross sections = 4.2387E-01 seconds\n",
" Total time in simulation = 1.4928E+00 seconds\n",
" Time in transport only = 1.3545E+00 seconds\n",
" Time in inactive batches = 1.3625E-01 seconds\n",
" Time in active batches = 1.3565E+00 seconds\n",
" Time synchronizing fission bank = 2.4053E-03 seconds\n",
" Sampling source sites = 1.6466E-03 seconds\n",
" SEND/RECV source sites = 5.6159E-04 seconds\n",
" Time accumulating tallies = 3.3647E-04 seconds\n",
" Total time for finalization = 1.6066E-02 seconds\n",
" Total time elapsed = 2.0336E+00 seconds\n",
" Calculation Rate (inactive) = 91743.2 neutrons/second\n",
" Calculation Rate (active) = 27644.5 neutrons/second\n",
" Total time for initialization = 3.4427e-01 seconds\n",
" Reading cross sections = 3.1628e-01 seconds\n",
" Total time in simulation = 3.7319e+00 seconds\n",
" Time in transport only = 3.6302e+00 seconds\n",
" Time in inactive batches = 4.9601e-01 seconds\n",
" Time in active batches = 3.2359e+00 seconds\n",
" Time synchronizing fission bank = 2.8100e-03 seconds\n",
" Sampling source sites = 2.4682e-03 seconds\n",
" SEND/RECV source sites = 3.2484e-04 seconds\n",
" Time accumulating tallies = 4.4538e-05 seconds\n",
" Total time for finalization = 9.3656e-04 seconds\n",
" Total time elapsed = 4.0859e+00 seconds\n",
" Calculation Rate (inactive) = 25201.2 particles/second\n",
" Calculation Rate (active) = 11588.7 particles/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.02889 +/- 0.00492\n",
" k-effective (Track-length) = 1.02842 +/- 0.00527\n",
" k-effective (Absorption) = 1.02637 +/- 0.00349\n",
" Combined k-effective = 1.02700 +/- 0.00291\n",
" Leakage Fraction = 0.01717 +/- 0.00107\n",
" k-effective (Collision) = 1.02889 +/- 0.00492\n",
" k-effective (Track-length) = 1.02842 +/- 0.00527\n",
" k-effective (Absorption) = 1.02637 +/- 0.00349\n",
" Combined k-effective = 1.02700 +/- 0.00291\n",
" Leakage Fraction = 0.01717 +/- 0.00107\n",
"\n"
]
}
@ -643,9 +574,8 @@
},
{
"cell_type": "code",
"execution_count": 22,
"execution_count": 20,
"metadata": {
"collapsed": true,
"scrolled": true
},
"outputs": [],
@ -665,13 +595,26 @@
},
{
"cell_type": "code",
"execution_count": 23,
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -699,7 +642,7 @@
"0 total (nu-fission / (absorption + current)) 1.02e+00 6.65e-03"
]
},
"execution_count": 23,
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
@ -729,13 +672,26 @@
},
{
"cell_type": "code",
"execution_count": 24,
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -770,7 +726,7 @@
"0 ((absorption + current) / (absorption + current)) 6.94e-01 4.61e-03 "
]
},
"execution_count": 24,
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
@ -794,13 +750,26 @@
},
{
"cell_type": "code",
"execution_count": 25,
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -835,7 +804,7 @@
"0 1.20e+00 9.61e-03 "
]
},
"execution_count": 25,
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
@ -858,13 +827,26 @@
},
{
"cell_type": "code",
"execution_count": 26,
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -901,7 +883,7 @@
"0 7.49e-01 6.09e-03 "
]
},
"execution_count": 26,
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
@ -922,13 +904,26 @@
},
{
"cell_type": "code",
"execution_count": 27,
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -965,7 +960,7 @@
"0 1.66e+00 1.44e-02 "
]
},
"execution_count": 27,
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
@ -985,13 +980,26 @@
},
{
"cell_type": "code",
"execution_count": 28,
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1026,7 +1034,7 @@
"0 ((absorption + current) / (absorption + current)) 9.85e-01 5.51e-03 "
]
},
"execution_count": 28,
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
@ -1045,13 +1053,26 @@
},
{
"cell_type": "code",
"execution_count": 29,
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1086,7 +1107,7 @@
"0 (absorption / (absorption + current)) 9.97e-01 7.55e-03 "
]
},
"execution_count": 29,
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
@ -1105,13 +1126,26 @@
},
{
"cell_type": "code",
"execution_count": 30,
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1148,7 +1182,7 @@
"0 (((((((absorption + current) / (absorption + c... 1.02e+00 1.88e-02 "
]
},
"execution_count": 30,
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
@ -1169,9 +1203,8 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": 29,
"metadata": {
"collapsed": true,
"scrolled": true
},
"outputs": [],
@ -1185,13 +1218,26 @@
},
{
"cell_type": "code",
"execution_count": 32,
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1312,7 +1358,7 @@
"7 (scatter / flux) 3.36e-03 1.34e-05 "
]
},
"execution_count": 32,
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
@ -1331,18 +1377,18 @@
},
{
"cell_type": "code",
"execution_count": 33,
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 6.65948580e-07]\n",
" [ 3.56632881e-01]]\n",
"[[[6.65948580e-07]\n",
" [3.56632881e-01]]\n",
"\n",
" [[ 7.25130446e-03]\n",
" [ 7.92016892e-03]]]\n"
" [[7.25130446e-03]\n",
" [7.92016892e-03]]]\n"
]
}
],
@ -1361,16 +1407,16 @@
},
{
"cell_type": "code",
"execution_count": 34,
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.00555547]]\n",
"[[[0.00555547]]\n",
"\n",
" [[ 0.00335828]]]\n"
" [[0.00335828]]]\n"
]
}
],
@ -1383,15 +1429,15 @@
},
{
"cell_type": "code",
"execution_count": 35,
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.22726611]\n",
" [ 0.00335828]]]\n"
"[[[0.22726611]\n",
" [0.00335828]]]\n"
]
}
],
@ -1412,13 +1458,26 @@
},
{
"cell_type": "code",
"execution_count": 36,
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1491,7 +1550,7 @@
"3 5.98e-04 "
]
},
"execution_count": 36,
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
@ -1504,13 +1563,26 @@
},
{
"cell_type": "code",
"execution_count": 37,
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1643,7 +1715,7 @@
"8 2.90e-03 "
]
},
"execution_count": 37,
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
@ -1673,7 +1745,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
"version": "3.7.0"
}
},
"nbformat": 4,

62
examples/jupyter/triso.ipynb
View file

File diff suppressed because one or more lines are too long

26
include/openmc/capi.h
View file

@ -10,15 +10,21 @@ extern "C" {
#endif
int openmc_calculate_volumes();
int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
int openmc_cell_filter_get_bins(int32_t index, const int32_t** cells, int32_t* n);
int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
int openmc_cell_get_id(int32_t index, int32_t* id);
int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T);
int openmc_cell_get_name(int32_t index, const char** name);
int openmc_cell_set_name(int32_t index, const char* name);
int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
int openmc_cell_set_id(int32_t index, int32_t id);
int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance);
int openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n);
int openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies);
int openmc_energy_filter_get_bins(int32_t index, const double** energies, size_t* n);
int openmc_energy_filter_set_bins(int32_t index, size_t n, const double* energies);
int openmc_energyfunc_filter_get_energy(int32_t index, size_t* n, const double** energy);
int openmc_energyfunc_filter_get_y(int32_t index, size_t* n, const double** y);
int openmc_energyfunc_filter_set_data(int32_t index, size_t n,
const double* energies, const double* y);
int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end);
@ -29,6 +35,8 @@ extern "C" {
int openmc_filter_set_id(int32_t index, int32_t id);
int openmc_finalize();
int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance);
int openmc_cell_bounding_box(const int32_t index, double* llc, double* urc);
int openmc_global_bounding_box(double* llc, double* urc);
int openmc_fission_bank(void** ptr, int64_t* n);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
@ -48,7 +56,7 @@ extern "C" {
int openmc_legendre_filter_set_order(int32_t index, int order);
int openmc_load_nuclide(const char* name);
int openmc_material_add_nuclide(int32_t index, const char name[], double density);
int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n);
int openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n);
int openmc_material_get_id(int32_t index, int32_t* id);
int openmc_material_get_fissionable(int32_t index, bool* fissionable);
int openmc_material_get_density(int32_t index, double* density);
@ -56,9 +64,11 @@ extern "C" {
int openmc_material_set_density(int32_t index, double density, const char* units);
int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density);
int openmc_material_set_id(int32_t index, int32_t id);
int openmc_material_get_name(int32_t index, const char** name);
int openmc_material_set_name(int32_t index, const char* name);
int openmc_material_set_volume(int32_t index, double volume);
int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n);
int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins);
int openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n);
int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins);
int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_mesh_get_id(int32_t index, int32_t* id);
@ -95,7 +105,7 @@ extern "C" {
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int* estimator);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n);
int openmc_tally_get_filters(int32_t index, const int32_t** indices, size_t* n);
int openmc_tally_get_n_realizations(int32_t index, int32_t* n);
int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
int openmc_tally_get_scores(int32_t index, int** scores, int* n);
@ -104,7 +114,7 @@ extern "C" {
int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]);
int openmc_tally_set_active(int32_t index, bool active);
int openmc_tally_set_estimator(int32_t index, const char* estimator);
int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices);
int openmc_tally_set_filters(int32_t index, size_t n, const int32_t* indices);
int openmc_tally_set_id(int32_t index, int32_t id);
int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
int openmc_tally_set_scores(int32_t index, int n, const char** scores);

121
include/openmc/cell.h
View file

@ -65,6 +65,8 @@ public:
//! \param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
BoundingBox bounding_box() const;
std::unique_ptr<UniversePartitioner> partitioner_;
};
@ -72,9 +74,77 @@ public:
//! A geometry primitive that links surfaces, universes, and materials
//==============================================================================
class Cell
{
class Cell {
public:
//----------------------------------------------------------------------------
// Constructors, destructors, factory functions
explicit Cell(pugi::xml_node cell_node);
Cell() {};
virtual ~Cell() = default;
//----------------------------------------------------------------------------
// Methods
//! \brief Determine if a cell contains the particle at a given location.
//!
//! The bounds of the cell are detemined by a logical expression involving
//! surface half-spaces. At initialization, the expression was converted
//! to RPN notation.
//!
//! The function is split into two cases, one for simple cells (those
//! involving only the intersection of half-spaces) and one for complex cells.
//! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
//! as we know that one half-space is not satisfied, we can exit. This
//! provides a performance benefit for the common case. In
//! contains_complex, we evaluate the RPN expression using a stack, similar to
//! how a RPN calculator would work.
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
virtual bool
contains(Position r, Direction u, int32_t on_surface) const = 0;
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
//! Get the BoundingBox for this cell.
virtual BoundingBox bounding_box() const = 0;
//----------------------------------------------------------------------------
// Accessors
//! Get the temperature of a cell instance
//! \param[in] instance Instance index. If -1 is given, the temperature for
//! the first instance is returned.
//! \return Temperature in [K]
double temperature(int32_t instance = -1) const;
//! Set the temperature of a cell instance
//! \param[in] T Temperature in [K]
//! \param[in] instance Instance index. If -1 is given, the temperature for
//! all instances is set.
void set_temperature(double T, int32_t instance = -1);
//! Get the name of a cell
//! \return Cell name
const std::string& name() const { return name_; };
//! Set the temperature of a cell instance
//! \param[in] name Cell name
void set_name(const std::string& name) { name_ = name; };
//----------------------------------------------------------------------------
// Data members
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
int type_; //!< Material, universe, or lattice
@ -116,41 +186,6 @@ public:
std::vector<double> rotation_;
std::vector<int32_t> offset_; //!< Distribcell offset table
explicit Cell(pugi::xml_node cell_node);
Cell() {};
//! \brief Determine if a cell contains the particle at a given location.
//!
//! The bounds of the cell are detemined by a logical expression involving
//! surface half-spaces. At initialization, the expression was converted
//! to RPN notation.
//!
//! The function is split into two cases, one for simple cells (those
//! involving only the intersection of half-spaces) and one for complex cells.
//! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
//! as we know that one half-space is not satisfied, we can exit. This
//! provides a performance benefit for the common case. In
//! contains_complex, we evaluate the RPN expression using a stack, similar to
//! how a RPN calculator would work.
//! \param r The 3D Cartesian coordinate to check.
//! \param u A direction used to "break ties" the coordinates are very
//! close to a surface.
//! \param on_surface The signed index of a surface that the coordinate is
//! known to be on. This index takes precedence over surface sense
//! calculations.
virtual bool
contains(Position r, Direction u, int32_t on_surface) const = 0;
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! @param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const = 0;
virtual ~Cell() {}
};
//==============================================================================
@ -170,9 +205,14 @@ public:
void to_hdf5(hid_t group_id) const;
BoundingBox bounding_box() const;
protected:
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
BoundingBox bounding_box_simple() const;
static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);
static void apply_demorgan(std::vector<int32_t>& rpn);
};
//==============================================================================
@ -181,16 +221,19 @@ protected:
class DAGCell : public Cell
{
public:
moab::DagMC* dagmc_ptr_;
DAGCell();
int32_t dag_index_;
bool contains(Position r, Direction u, int32_t on_surface) const;
std::pair<double, int32_t>
distance(Position r, Direction u, int32_t on_surface) const;
BoundingBox bounding_box() const;
void to_hdf5(hid_t group_id) const;
moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of cell
};
#endif

81
include/openmc/material.h
View file

@ -6,6 +6,7 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include <hdf5.h>
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
@ -35,6 +36,7 @@ extern std::unordered_map<int32_t, int32_t> material_map;
class Material
{
public:
//----------------------------------------------------------------------------
// Types
struct ThermalTable {
int index_table; //!< Index of table in data::thermal_scatt
@ -42,11 +44,15 @@ public:
double fraction; //!< How often to use table
};
// Constructors
//----------------------------------------------------------------------------
// Constructors, destructors, factory functions
Material() {};
explicit Material(pugi::xml_node material_node);
~Material();
//----------------------------------------------------------------------------
// Methods
void calculate_xs(Particle& p) const;
//! Assign thermal scattering tables to specific nuclides within the material
@ -59,12 +65,72 @@ public:
//! Finalize the material, assigning tables, normalize density, etc.
void finalize();
//! Set total density of the material
int set_density(double density, std::string units);
//! Write material data to HDF5
void to_hdf5(hid_t group) const;
//! Add nuclide to the material
//
//! \param[in] nuclide Name of the nuclide
//! \param[in] density Density of the nuclide in [atom/b-cm]
void add_nuclide(const std::string& nuclide, double density);
//! Set atom densities for the material
//
//! \param[in] name Name of each nuclide
//! \param[in] density Density of each nuclide in [atom/b-cm]
void set_densities(const std::vector<std::string>& name,
const std::vector<double>& density);
//----------------------------------------------------------------------------
// Accessors
//! Get density in [atom/b-cm]
//! \return Density in [atom/b-cm]
double density() const { return density_; }
//! Get density in [g/cm^3]
//! \return Density in [g/cm^3]
double density_gpcc() const { return density_gpcc_; }
//! Get name
//! \return Material name
const std::string& name() const { return name_; }
//! Set name
void set_name(const std::string& name) { name_ = name; }
//! Set total density of the material
//
//! \param[in] density Density value
//! \param[in] units Units of density
void set_density(double density, gsl::cstring_span units);
//! Get nuclides in material
//! \return Indices into the global nuclides vector
gsl::span<const int> nuclides() const { return {nuclide_.data(), nuclide_.size()}; }
//! Get densities of each nuclide in material
//! \return Densities in [atom/b-cm]
gsl::span<const double> densities() const { return {atom_density_.data(), atom_density_.size()}; }
//! Get ID of material
//! \return ID of material
int32_t id() const { return id_; }
//! Assign a unique ID to the material
//! \param[in] Unique ID to assign. A value of -1 indicates that an ID
//! should be automatically assigned.
void set_id(int32_t id);
//! Get whether material is fissionable
//! \return Whether material is fissionable
bool fissionable() const { return fissionable_; }
//! Get volume of material
//! \return Volume in [cm^3]
double volume() const;
//----------------------------------------------------------------------------
// Data
int32_t id_; //!< Unique ID
std::string name_; //!< Name of material
@ -95,6 +161,9 @@ public:
std::unique_ptr<Bremsstrahlung> ttb_;
private:
//----------------------------------------------------------------------------
// Private methods
//! Calculate the collision stopping power
void collision_stopping_power(double* s_col, bool positron);
@ -106,6 +175,10 @@ private:
void calculate_neutron_xs(Particle& p) const;
void calculate_photon_xs(Particle& p) const;
//----------------------------------------------------------------------------
// Private data members
gsl::index index_;
};
//==============================================================================

68
include/openmc/surface.h
View file

@ -43,12 +43,46 @@ namespace model {
struct BoundingBox
{
double xmin;
double xmax;
double ymin;
double ymax;
double zmin;
double zmax;
double xmin = -INFTY;
double xmax = INFTY;
double ymin = -INFTY;
double ymax = INFTY;
double zmin = -INFTY;
double zmax = INFTY;
inline BoundingBox operator &(const BoundingBox& other) {
BoundingBox result = *this;
return result &= other;
}
inline BoundingBox operator |(const BoundingBox& other) {
BoundingBox result = *this;
return result |= other;
}
// intersect operator
inline BoundingBox& operator &=(const BoundingBox& other) {
xmin = std::max(xmin, other.xmin);
xmax = std::min(xmax, other.xmax);
ymin = std::max(ymin, other.ymin);
ymax = std::min(ymax, other.ymax);
zmin = std::max(zmin, other.zmin);
zmax = std::min(zmax, other.zmax);
return *this;
}
// union operator
inline BoundingBox& operator |=(const BoundingBox& other) {
xmin = std::min(xmin, other.xmin);
xmax = std::max(xmax, other.xmax);
ymin = std::min(ymin, other.ymin);
ymax = std::max(ymax, other.ymax);
zmin = std::min(zmin, other.zmin);
zmax = std::max(zmax, other.zmax);
return *this;
}
};
//==============================================================================
@ -105,6 +139,8 @@ public:
//TODO: this probably needs to include i_periodic for PeriodicSurface
virtual void to_hdf5(hid_t group_id) const = 0;
//! Get the BoundingBox for this surface.
virtual BoundingBox bounding_box(bool pos_side) const { return {}; }
};
class CSGSurface : public Surface
@ -126,18 +162,17 @@ protected:
class DAGSurface : public Surface
{
public:
moab::DagMC* dagmc_ptr_;
DAGSurface();
int32_t dag_index_;
double evaluate(Position r) const;
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
Direction reflect(Position r, Direction u) const;
//! Get the bounding box of this surface.
BoundingBox bounding_box() const;
void to_hdf5(hid_t group_id) const;
moab::DagMC* dagmc_ptr_; //!< Pointer to DagMC instance
int32_t dag_index_; //!< DagMC index of surface
};
#endif
//==============================================================================
@ -166,8 +201,6 @@ public:
virtual bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const = 0;
//! Get the bounding box for this surface.
virtual BoundingBox bounding_box() const = 0;
};
//==============================================================================
@ -186,7 +219,7 @@ public:
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
BoundingBox bounding_box(bool pos_side) const;
double x0_;
};
@ -207,7 +240,7 @@ public:
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
BoundingBox bounding_box(bool pos_side) const;
double y0_;
};
@ -228,7 +261,7 @@ public:
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
BoundingBox bounding_box(bool pos_side) const;
double z0_;
};
@ -249,7 +282,6 @@ public:
void to_hdf5_inner(hid_t group_id) const;
bool periodic_translate(const PeriodicSurface* other, Position& r,
Direction& u) const;
BoundingBox bounding_box() const;
double A_, B_, C_, D_;
};
@ -269,6 +301,7 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
BoundingBox bounding_box(bool pos_side) const;
double y0_, z0_, radius_;
};
@ -288,6 +321,7 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
BoundingBox bounding_box(bool pos_side) const;
double x0_, z0_, radius_;
};
@ -307,6 +341,7 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
BoundingBox bounding_box(bool pos_side) const;
double x0_, y0_, radius_;
};
@ -326,6 +361,7 @@ public:
double distance(Position r, Direction u, bool coincident) const;
Direction normal(Position r) const;
void to_hdf5_inner(hid_t group_id) const;
BoundingBox bounding_box(bool pos_side) const;
double x0_, y0_, z0_, radius_;
};

54
include/openmc/tallies/filter.h
View file

@ -7,6 +7,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
#include "pugixml.hpp"
@ -43,13 +45,34 @@ namespace openmc {
class Filter
{
public:
virtual ~Filter() = default;
//----------------------------------------------------------------------------
// Constructors, destructors, factory functions
virtual std::string type() const = 0;
Filter();
virtual ~Filter();
//! Create a new tally filter
//
//! \param[in] type Type of the filter
//! \param[in] id Unique ID for the filter. If none is passed, an ID is
//! automatically assigned
//! \return Pointer to the new filter object
static Filter* create(const std::string& type, int32_t id = -1);
//! Create a new tally filter from an XML node
//
//! \param[in] node XML node
//! \return Pointer to the new filter object
static Filter* create(pugi::xml_node node);
//! Uses an XML input to fill the filter's data fields.
virtual void from_xml(pugi::xml_node node) = 0;
//----------------------------------------------------------------------------
// Methods
virtual std::string type() const = 0;
//! Matches a tally event to a set of filter bins and weights.
//!
//! \param[out] match will contain the matching bins and corresponding
@ -71,11 +94,32 @@ public:
//! "Incoming Energy [0.625E-6, 20.0)".
virtual std::string text_label(int bin) const = 0;
virtual void initialize() {}
//----------------------------------------------------------------------------
// Accessors
int32_t id_;
//! Get unique ID of filter
//! \return Unique ID
int32_t id() const { return id_; }
//! Assign a unique ID to the filter
//! \param[in] Unique ID to assign. A value of -1 indicates that an ID should
//! be automatically assigned
void set_id(int32_t id);
//! Get number of bins
//! \return Number of bins
int n_bins() const { return n_bins_; }
gsl::index index() const { return index_; }
//----------------------------------------------------------------------------
// Data members
protected:
int n_bins_;
private:
int32_t id_ {-1};
gsl::index index_;
};
//==============================================================================
@ -99,8 +143,6 @@ namespace model {
// Non-member functions
//==============================================================================
Filter* allocate_filter(const std::string& type);
//! Make sure index corresponds to a valid filter
int verify_filter(int32_t index);

17
include/openmc/tallies/filter_azimuthal.h
View file

@ -4,6 +4,8 @@
#include <string>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -15,8 +17,14 @@ namespace openmc {
class AzimuthalFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~AzimuthalFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "azimuthal";}
void from_xml(pugi::xml_node node) override;
@ -28,6 +36,15 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_bins(gsl::span<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

21
include/openmc/tallies/filter_cell.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -16,14 +18,18 @@ namespace openmc {
class CellFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~CellFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "cell";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
@ -31,6 +37,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const std::vector<int32_t>& cells() const { return cells_; }
void set_cells(gsl::span<int32_t> cells);
protected:
//----------------------------------------------------------------------------
// Data members
//! The indices of the cells binned by this filter.
std::vector<int32_t> cells_;

3
include/openmc/tallies/filter_cellborn.h
View file

@ -14,6 +14,9 @@ namespace openmc {
class CellbornFilter : public CellFilter
{
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "cellborn";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)

3
include/openmc/tallies/filter_cellfrom.h
View file

@ -14,6 +14,9 @@ namespace openmc {
class CellFromFilter : public CellFilter
{
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "cellfrom";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)

19
include/openmc/tallies/filter_delayedgroup.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -17,8 +19,14 @@ namespace openmc {
class DelayedGroupFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~DelayedGroupFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "delayedgroup";}
void from_xml(pugi::xml_node node) override;
@ -30,6 +38,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const std::vector<int>& groups() const { return groups_; }
void set_groups(gsl::span<int> groups);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<int> groups_;
};

19
include/openmc/tallies/filter_distribcell.h
View file

@ -14,14 +14,18 @@ namespace openmc {
class DistribcellFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~DistribcellFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "distribcell";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
@ -29,6 +33,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
int32_t cell() const { return cell_; }
void set_cell(int32_t cell);
private:
//----------------------------------------------------------------------------
// Data members
int32_t cell_;
};

23
include/openmc/tallies/filter_energy.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -14,8 +16,14 @@ namespace openmc {
class EnergyFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~EnergyFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "energy";}
void from_xml(pugi::xml_node node) override;
@ -27,6 +35,18 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const std::vector<double>& bins() const { return bins_; }
void set_bins(gsl::span<const double> bins);
bool matches_transport_groups() const { return matches_transport_groups_; }
protected:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
//! True if transport group number can be used directly to get bin number
@ -43,6 +63,9 @@ public:
class EnergyoutFilter : public EnergyFilter
{
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "energyout";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)

17
include/openmc/tallies/filter_energyfunc.h
View file

@ -15,6 +15,9 @@ namespace openmc {
class EnergyFunctionFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
EnergyFunctionFilter()
: Filter {}
{
@ -23,6 +26,9 @@ public:
~EnergyFunctionFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "energyfunction";}
void from_xml(pugi::xml_node node) override;
@ -34,6 +40,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const std::vector<double>& energy() const { return energy_; }
const std::vector<double>& y() const { return y_; }
void set_data(gsl::span<const double> energy, gsl::span<const double> y);
private:
//----------------------------------------------------------------------------
// Data members
//! Incident neutron energy interpolation grid.
std::vector<double> energy_;

17
include/openmc/tallies/filter_legendre.h
View file

@ -14,8 +14,14 @@ namespace openmc {
class LegendreFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~LegendreFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "legendre";}
void from_xml(pugi::xml_node node) override;
@ -27,6 +33,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
int order() const { return order_; }
void set_order(int order);
private:
//----------------------------------------------------------------------------
// Data members
int order_;
};

23
include/openmc/tallies/filter_material.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -16,14 +18,18 @@ namespace openmc {
class MaterialFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~MaterialFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "material";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
@ -31,6 +37,19 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
std::vector<int32_t>& materials() { return materials_; }
const std::vector<int32_t>& materials() const { return materials_; }
void set_materials(gsl::span<const int32_t> materials);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the materials binned by this filter.
std::vector<int32_t> materials_;

12
include/openmc/tallies/filter_mesh.h
View file

@ -16,8 +16,14 @@ namespace openmc {
class MeshFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~MeshFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "mesh";}
void from_xml(pugi::xml_node node) override;
@ -29,11 +35,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
virtual int32_t mesh() const {return mesh_;}
virtual void set_mesh(int32_t mesh);
protected:
//----------------------------------------------------------------------------
// Data members
int32_t mesh_;
};

6
include/openmc/tallies/filter_meshsurface.h
View file

@ -8,6 +8,9 @@ namespace openmc {
class MeshSurfaceFilter : public MeshFilter
{
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "meshsurface";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
@ -15,6 +18,9 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_mesh(int32_t mesh) override;
};

17
include/openmc/tallies/filter_mu.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -15,8 +17,14 @@ namespace openmc {
class MuFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~MuFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "mu";}
void from_xml(pugi::xml_node node) override;
@ -28,6 +36,15 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_bins(gsl::span<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

17
include/openmc/tallies/filter_particle.h
View file

@ -15,8 +15,14 @@ namespace openmc {
class ParticleFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~ParticleFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "particle";}
void from_xml(pugi::xml_node node) override;
@ -28,6 +34,17 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const std::vector<Particle::Type>& particles() const { return particles_; }
void set_particles(gsl::span<Particle::Type> particles);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<Particle::Type> particles_;
};

17
include/openmc/tallies/filter_polar.h
View file

@ -4,6 +4,8 @@
#include <cmath>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -15,8 +17,14 @@ namespace openmc {
class PolarFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~PolarFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "polar";}
void from_xml(pugi::xml_node node) override;
@ -28,6 +36,15 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_bins(gsl::span<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

23
include/openmc/tallies/filter_sph_harm.h
View file

@ -3,6 +3,8 @@
#include <string>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -18,8 +20,14 @@ enum class SphericalHarmonicsCosine {
class SphericalHarmonicsFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~SphericalHarmonicsFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "sphericalharmonics";}
void from_xml(pugi::xml_node node) override;
@ -31,6 +39,21 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
int order() const { return order_; }
void set_order(int order);
SphericalHarmonicsCosine cosine() const { return cosine_; }
void set_cosine(gsl::cstring_span cosine);
private:
//----------------------------------------------------------------------------
// Data members
int order_;
//! The type of angle that this filter measures when binning events.

23
include/openmc/tallies/filter_sptl_legendre.h
View file

@ -18,8 +18,14 @@ enum class LegendreAxis {
class SpatialLegendreFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~SpatialLegendreFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "spatiallegendre";}
void from_xml(pugi::xml_node node) override;
@ -31,6 +37,23 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
int order() const { return order_; }
void set_order(int order);
LegendreAxis axis() const { return axis_; }
void set_axis(LegendreAxis axis);
double min() const { return min_; }
double max() const { return max_; }
void set_minmax(double min, double max);
private:
//----------------------------------------------------------------------------
// Data members
int order_;
//! The Cartesian coordinate axis that the Legendre expansion is applied to.

19
include/openmc/tallies/filter_surface.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -16,14 +18,18 @@ namespace openmc {
class SurfaceFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~SurfaceFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "surface";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
@ -31,6 +37,15 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_surfaces(gsl::span<int32_t> surfaces);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the surfaces binned by this filter.
std::vector<int32_t> surfaces_;

19
include/openmc/tallies/filter_universe.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include "openmc/tallies/filter.h"
namespace openmc {
@ -16,14 +18,18 @@ namespace openmc {
class UniverseFilter : public Filter
{
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~UniverseFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "universe";}
void from_xml(pugi::xml_node node) override;
void initialize() override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
@ -31,6 +37,15 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_universes(gsl::span<int32_t> universes);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the universes binned by this filter.
std::vector<int32_t> universes_;

32
include/openmc/tallies/filter_zernike.h
View file

@ -14,10 +14,16 @@ namespace openmc {
class ZernikeFilter : public Filter
{
public:
std::string type() const override {return "zernike";}
//----------------------------------------------------------------------------
// Constructors, destructors
~ZernikeFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "zernike";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
@ -27,10 +33,25 @@ public:
std::string text_label(int bin) const override;
int order() const {return order_;}
//----------------------------------------------------------------------------
// Accessors
int order() const { return order_; }
virtual void set_order(int order);
double x() const { return x_; }
void set_x(double x) { x_ = x; }
double y() const { return y_; }
void set_y(double y) { y_ = y; }
double r() const { return r_; }
void set_r(double r) { r_ = r; }
//----------------------------------------------------------------------------
// Data members
protected:
//! Cartesian x coordinate for the origin of this expansion.
double x_;
@ -40,7 +61,6 @@ public:
//! Maximum radius from the origin covered by this expansion.
double r_;
protected:
int order_;
};
@ -51,6 +71,9 @@ protected:
class ZernikeRadialFilter : public ZernikeFilter
{
public:
//----------------------------------------------------------------------------
// Methods
std::string type() const override {return "zernikeradial";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
@ -58,6 +81,9 @@ public:
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
void set_order(int order) override;
};

25
include/openmc/tallies/tally.h
View file

@ -2,8 +2,10 @@
#define OPENMC_TALLIES_TALLY_H
#include "openmc/constants.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/trigger.h"
#include <gsl/gsl>
#include "pugixml.hpp"
#include "xtensor/xfixed.hpp"
#include "xtensor/xtensor.hpp"
@ -21,24 +23,33 @@ namespace openmc {
class Tally {
public:
Tally();
//----------------------------------------------------------------------------
// Constructors, destructors, factory functions
explicit Tally(int32_t id);
explicit Tally(pugi::xml_node node);
~Tally();
static Tally* create(int32_t id = -1);
void init_from_xml(pugi::xml_node node);
//----------------------------------------------------------------------------
// Accessors
void set_id(int32_t id);
void set_active(bool active) { active_ = active; }
void set_scores(pugi::xml_node node);
void set_scores(std::vector<std::string> scores);
void set_scores(const std::vector<std::string>& scores);
void set_nuclides(pugi::xml_node node);
//----------------------------------------------------------------------------
// Methods for getting and setting filter/stride data.
void set_nuclides(const std::vector<std::string>& nuclides);
const std::vector<int32_t>& filters() const {return filters_;}
int32_t filters(int i) const {return filters_[i];}
void set_filters(const int32_t filter_indices[], int n);
void set_filters(gsl::span<Filter*> filters);
int32_t strides(int i) const {return strides_[i];}
@ -112,6 +123,8 @@ private:
std::vector<int32_t> strides_;
int32_t n_filter_bins_ {0};
gsl::index index_;
};
//==============================================================================

38
openmc/capi/cell.py
View file

@ -1,3 +1,5 @@
import sys
from collections.abc import Mapping, Iterable
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from weakref import WeakValueDictionary
@ -28,6 +30,12 @@ _dll.openmc_cell_get_temperature.argtypes = [
c_int32, POINTER(c_int32), POINTER(c_double)]
_dll.openmc_cell_get_temperature.restype = c_int
_dll.openmc_cell_get_temperature.errcheck = _error_handler
_dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)]
_dll.openmc_cell_get_name.restype = c_int
_dll.openmc_cell_get_name.errcheck = _error_handler
_dll.openmc_cell_set_name.argtypes = [c_int32, c_char_p]
_dll.openmc_cell_set_name.restype = c_int
_dll.openmc_cell_set_name.errcheck = _error_handler
_dll.openmc_cell_set_fill.argtypes = [
c_int32, c_int, c_int32, POINTER(c_int32)]
_dll.openmc_cell_set_fill.restype = c_int
@ -43,6 +51,11 @@ _dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_cell_index.restype = c_int
_dll.openmc_get_cell_index.errcheck = _error_handler
_dll.cells_size.restype = c_int
_dll.openmc_cell_bounding_box.argtypes = [c_int,
POINTER(c_double),
POINTER(c_double)]
_dll.openmc_cell_bounding_box.restype = c_int
_dll.openmc_cell_bounding_box.errcheck = _error_handler
class Cell(_FortranObjectWithID):
@ -102,6 +115,17 @@ class Cell(_FortranObjectWithID):
def id(self, cell_id):
_dll.openmc_cell_set_id(self._index, cell_id)
@property
def name(self):
name = c_char_p()
_dll.openmc_cell_get_name(self._index, name)
return name.value.decode()
@name.setter
def name(self, name):
name_ptr = c_char_p(name.encode())
_dll.openmc_cell_set_name(self._index, name_ptr)
@property
def fill(self):
fill_type = c_int()
@ -166,6 +190,20 @@ class Cell(_FortranObjectWithID):
_dll.openmc_cell_set_temperature(self._index, T, instance)
@property
def bounding_box(self):
inf = sys.float_info.max
llc = np.zeros(3)
urc = np.zeros(3)
_dll.openmc_cell_bounding_box(self._index,
llc.ctypes.data_as(POINTER(c_double)),
urc.ctypes.data_as(POINTER(c_double)))
llc[llc == inf] = np.inf
urc[urc == inf] = np.inf
llc[llc == -inf] = -np.inf
urc[urc == -inf] = -np.inf
return llc, urc
class _CellMapping(Mapping):
def __getitem__(self, key):

20
openmc/capi/core.py
View file

@ -1,3 +1,5 @@
import sys
from contextlib import contextmanager
from ctypes import (CDLL, c_bool, c_int, c_int32, c_int64, c_double, c_char_p,
c_char, POINTER, Structure, c_void_p, create_string_buffer)
@ -73,8 +75,26 @@ _dll.openmc_simulation_finalize.errcheck = _error_handler
_dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)]
_dll.openmc_statepoint_write.restype = c_int
_dll.openmc_statepoint_write.errcheck = _error_handler
_dll.openmc_global_bounding_box.argtypes = [POINTER(c_double),
POINTER(c_double)]
_dll.openmc_global_bounding_box.restype = c_int
_dll.openmc_global_bounding_box.errcheck = _error_handler
def global_bounding_box():
"""Calculate a global bounding box for the model"""
inf = sys.float_info.max
llc = np.zeros(3)
urc = np.zeros(3)
_dll.openmc_global_bounding_box(llc.ctypes.data_as(POINTER(c_double)),
urc.ctypes.data_as(POINTER(c_double)))
llc[llc == inf] = np.inf
urc[urc == inf] = np.inf
llc[llc == -inf] = -np.inf
urc[urc == -inf] = -np.inf
return llc, urc
def calculate_volumes():
"""Run stochastic volume calculation"""
_dll.openmc_calculate_volumes()

4
openmc/capi/error.py
View file

@ -36,4 +36,6 @@ def _error_handler(err, func, args):
elif err == errcode('OPENMC_E_WARNING'):
warn(msg)
elif err < 0:
raise exc.OpenMCError("Unknown error encountered (code {}).".format(err))
if not msg:
msg = "Unknown error encountered (code {}).".format(err)
raise exc.OpenMCError(msg)

67
openmc/capi/filter.py
View file

@ -28,12 +28,24 @@ _dll.openmc_cell_filter_get_bins.argtypes = [
_dll.openmc_cell_filter_get_bins.restype = c_int
_dll.openmc_cell_filter_get_bins.errcheck = _error_handler
_dll.openmc_energy_filter_get_bins.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(c_int32)]
c_int32, POINTER(POINTER(c_double)), POINTER(c_size_t)]
_dll.openmc_energy_filter_get_bins.restype = c_int
_dll.openmc_energy_filter_get_bins.errcheck = _error_handler
_dll.openmc_energy_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_double)]
_dll.openmc_energy_filter_set_bins.argtypes = [c_int32, c_size_t, POINTER(c_double)]
_dll.openmc_energy_filter_set_bins.restype = c_int
_dll.openmc_energy_filter_set_bins.errcheck = _error_handler
_dll.openmc_energyfunc_filter_set_data.restype = c_int
_dll.openmc_energyfunc_filter_set_data.errcheck = _error_handler
_dll.openmc_energyfunc_filter_set_data.argtypes = [
c_int32, c_size_t, POINTER(c_double), POINTER(c_double)]
_dll.openmc_energyfunc_filter_get_energy.resttpe = c_int
_dll.openmc_energyfunc_filter_get_energy.errcheck = _error_handler
_dll.openmc_energyfunc_filter_get_energy.argtypes = [
c_int32, POINTER(c_size_t), POINTER(POINTER(c_double))]
_dll.openmc_energyfunc_filter_get_y.resttpe = c_int
_dll.openmc_energyfunc_filter_get_y.errcheck = _error_handler
_dll.openmc_energyfunc_filter_get_y.argtypes = [
c_int32, POINTER(c_size_t), POINTER(POINTER(c_double))]
_dll.openmc_filter_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_filter_get_id.restype = c_int
_dll.openmc_filter_get_id.errcheck = _error_handler
@ -53,10 +65,10 @@ _dll.openmc_legendre_filter_set_order.argtypes = [c_int32, c_int]
_dll.openmc_legendre_filter_set_order.restype = c_int
_dll.openmc_legendre_filter_set_order.errcheck = _error_handler
_dll.openmc_material_filter_get_bins.argtypes = [
c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)]
c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)]
_dll.openmc_material_filter_get_bins.restype = c_int
_dll.openmc_material_filter_get_bins.errcheck = _error_handler
_dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_int32)]
_dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_size_t, POINTER(c_int32)]
_dll.openmc_material_filter_set_bins.restype = c_int
_dll.openmc_material_filter_set_bins.errcheck = _error_handler
_dll.openmc_mesh_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
@ -94,6 +106,7 @@ _dll.openmc_zernike_filter_set_order.restype = c_int
_dll.openmc_zernike_filter_set_order.errcheck = _error_handler
_dll.tally_filters_size.restype = c_size_t
class Filter(_FortranObjectWithID):
__instances = WeakValueDictionary()
@ -148,7 +161,7 @@ class EnergyFilter(Filter):
@property
def bins(self):
energies = POINTER(c_double)()
n = c_int32()
n = c_size_t()
_dll.openmc_energy_filter_get_bins(self._index, energies, n)
return as_array(energies, (n.value,))
@ -200,6 +213,48 @@ class DistribcellFilter(Filter):
class EnergyFunctionFilter(Filter):
filter_type = 'energyfunction'
def __new__(cls, energy=None, y=None, uid=None, new=True, index=None):
return super().__new__(cls, uid=uid, new=new, index=index)
def __init__(self, energy=None, y=None, uid=None, new=True, index=None):
if (energy is None) != (y is None):
raise AttributeError("Need both energy and y or neither")
super().__init__(uid, new, index)
if energy is not None:
self.set_data(energy, y)
def set_data(self, energy, y):
"""Set the interpolation information for the filter
Parameters
----------
energy : numpy.ndarray
Independent variable for the interpolation
y : numpy.ndarray
Dependent variable for the interpolation
"""
energy_array = np.asarray(energy)
y_array = np.asarray(y)
energy_p = energy_array.ctypes.data_as(POINTER(c_double))
y_p = y_array.ctypes.data_as(POINTER(c_double))
_dll.openmc_energyfunc_filter_set_data(
self._index, len(energy_array), energy_p, y_p)
@property
def energy(self):
return self._get_attr(_dll.openmc_energyfunc_filter_get_energy)
@property
def y(self):
return self._get_attr(_dll.openmc_energyfunc_filter_get_y)
def _get_attr(self, cfunc):
array_p = POINTER(c_double)()
n = c_size_t()
cfunc(self._index, n, array_p)
return as_array(array_p, (n.value, ))
class LegendreFilter(Filter):
filter_type = 'legendre'
@ -231,7 +286,7 @@ class MaterialFilter(Filter):
@property
def bins(self):
materials = POINTER(c_int32)()
n = c_int32()
n = c_size_t()
_dll.openmc_material_filter_get_bins(self._index, materials, n)
return [Material(index=materials[i]) for i in range(n.value)]

17
openmc/capi/material.py
View file

@ -48,6 +48,12 @@ _dll.openmc_material_set_densities.errcheck = _error_handler
_dll.openmc_material_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_material_set_id.restype = c_int
_dll.openmc_material_set_id.errcheck = _error_handler
_dll.openmc_material_get_name.argtypes = [c_int32, POINTER(c_char_p)]
_dll.openmc_material_get_name.restype = c_int
_dll.openmc_material_get_name.errcheck = _error_handler
_dll.openmc_material_set_name.argtypes = [c_int32, c_char_p]
_dll.openmc_material_set_name.restype = c_int
_dll.openmc_material_set_name.errcheck = _error_handler
_dll.openmc_material_set_volume.argtypes = [c_int32, c_double]
_dll.openmc_material_set_volume.restype = c_int
_dll.openmc_material_set_volume.errcheck = _error_handler
@ -124,6 +130,17 @@ class Material(_FortranObjectWithID):
def id(self, mat_id):
_dll.openmc_material_set_id(self._index, mat_id)
@property
def name(self):
name = c_char_p()
_dll.openmc_material_get_name(self._index, name)
return name.value.decode()
@name.setter
def name(self, name):
name_ptr = c_char_p(name.encode())
_dll.openmc_material_set_name(self._index, name_ptr)
@property
def volume(self):
volume = c_double()

6
openmc/capi/tally.py
View file

@ -36,7 +36,7 @@ _dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_id.restype = c_int
_dll.openmc_tally_get_id.errcheck = _error_handler
_dll.openmc_tally_get_filters.argtypes = [
c_int32, POINTER(POINTER(c_int32)), POINTER(c_int)]
c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)]
_dll.openmc_tally_get_filters.restype = c_int
_dll.openmc_tally_get_filters.errcheck = _error_handler
_dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)]
@ -63,7 +63,7 @@ _dll.openmc_tally_results.errcheck = _error_handler
_dll.openmc_tally_set_active.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_active.restype = c_int
_dll.openmc_tally_set_active.errcheck = _error_handler
_dll.openmc_tally_set_filters.argtypes = [c_int32, c_int, POINTER(c_int32)]
_dll.openmc_tally_set_filters.argtypes = [c_int32, c_size_t, POINTER(c_int32)]
_dll.openmc_tally_set_filters.restype = c_int
_dll.openmc_tally_set_filters.errcheck = _error_handler
_dll.openmc_tally_set_estimator.argtypes = [c_int32, c_char_p]
@ -249,7 +249,7 @@ class Tally(_FortranObjectWithID):
@property
def filters(self):
filt_idx = POINTER(c_int32)()
n = c_int()
n = c_size_t()
_dll.openmc_tally_get_filters(self._index, filt_idx, n)
return [_get_filter(filt_idx[i]) for i in range(n.value)]

2
openmc/data/ace.py
View file

@ -411,7 +411,7 @@ class Library(EqualityMixin):
# after it). If it's too short, then we apply the ENDF float regular
# expression. We don't do this by default because it's expensive!
if xss.size != nxs[1] + 1:
datastr = ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
datastr = ENDF_FLOAT_RE.sub(r'\1e\2\3', datastr)
xss = np.fromstring(datastr, sep=' ')
assert xss.size == nxs[1] + 1

13
openmc/data/neutron.py
View file

@ -801,15 +801,14 @@ class IncidentNeutron(EqualityMixin):
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
pendf_file = os.path.join(tmpdir, 'pendf')
kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
kwargs.setdefault('pendf', os.path.join(tmpdir, 'pendf'))
kwargs['evaluation'] = evaluation
make_ace(filename, temperatures, ace_file, xsdir_file,
pendf_file, **kwargs)
make_ace(filename, temperatures, **kwargs)
# Create instance from ACE tables within library
lib = Library(ace_file)
lib = Library(kwargs['ace'])
data = cls.from_ace(lib.tables[0])
for table in lib.tables[1:]:
data.add_temperature_from_ace(table)
@ -821,7 +820,7 @@ class IncidentNeutron(EqualityMixin):
# Add 0K elastic scattering cross section
if '0K' not in data.energy:
pendf = Evaluation(pendf_file)
pendf = Evaluation(kwargs['pendf'])
file_obj = StringIO(pendf.section[3, 2])
get_head_record(file_obj)
params, xs = get_tab1_record(file_obj)

9
openmc/data/thermal.py
View file

@ -761,15 +761,14 @@ class ThermalScattering(EqualityMixin):
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
kwargs['evaluation'] = evaluation
kwargs['evaluation_thermal'] = evaluation_thermal
make_ace_thermal(filename, filename_thermal, temperatures,
ace_file, xsdir_file, **kwargs)
make_ace_thermal(filename, filename_thermal, temperatures, **kwargs)
# Create instance from ACE tables within library
lib = Library(ace_file)
lib = Library(kwargs['ace'])
data = cls.from_ace(lib.tables[0])
for table in lib.tables[1:]:
data.add_temperature_from_ace(table)

195
openmc/deplete/abc.py
View file

@ -7,11 +7,15 @@ to run a full depletion simulation.
from collections import namedtuple
import os
from pathlib import Path
from abc import ABCMeta, abstractmethod
from abc import ABC, abstractmethod
from xml.etree import ElementTree as ET
from warnings import warn
from numbers import Real
from openmc.data import DataLibrary
from numpy import nonzero, empty
from openmc.data import DataLibrary, JOULE_PER_EV
from openmc.checkvalue import check_type, check_greater_than
from .chain import Chain
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
@ -34,7 +38,7 @@ except AttributeError:
pass
class TransportOperator(metaclass=ABCMeta):
class TransportOperator(ABC):
"""Abstract class defining a transport operator
Each depletion integrator is written to work with a generic transport
@ -52,17 +56,21 @@ class TransportOperator(metaclass=ABCMeta):
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``.
dilute_initial : float, optional
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
Defaults to 1.0e3.
Attributes
----------
dilute_initial : float
Initial atom density to add for nuclides that are zero in initial
condition to ensure they exist in the decay chain. Only done for
nuclides with reaction rates. Defaults to 1.0e3.
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
"""
def __init__(self, chain_file=None, fission_q=None):
self.dilute_initial = 1.0e3
def __init__(self, chain_file=None, fission_q=None, dilute_initial=1.0e3):
self.dilute_initial = dilute_initial
self.output_dir = '.'
# Read depletion chain
@ -86,6 +94,17 @@ class TransportOperator(metaclass=ABCMeta):
FutureWarning)
self.chain = Chain.from_xml(chain_file, fission_q)
@property
def dilute_initial(self):
"""Initial atom density for nuclides with zero initial concentration"""
return self._dilute_initial
@dilute_initial.setter
def dilute_initial(self, value):
check_type("dilute_initial", value, Real)
check_greater_than("dilute_initial", value, 0.0, equality=True)
self._dilute_initial = value
@abstractmethod
def __call__(self, vec):
"""Runs a simulation.
@ -158,3 +177,161 @@ class TransportOperator(metaclass=ABCMeta):
def finalize(self):
pass
class ReactionRateHelper(ABC):
"""Abstract class for generating reaction rates for operators
Responsible for generating reaction rate tallies for burnable
materials, given nuclides and scores from the operator.
Reaction rates are passed back to the operator for be used in
an :class:`openmc.deplete.OperatorResult` instance
Parameters
----------
n_nucs : int
Number of burnable nuclides tracked by :class:`openmc.deplete.Operator`
n_react : int
Number of reactions tracked by :class:`openmc.deplete.Operator`
Attributes
----------
nuclides : list of str
All nuclides with desired reaction rates.
"""
def __init__(self, n_nucs, n_react):
self._nuclides = None
self._rate_tally = None
self._results_cache = empty((n_nucs, n_react))
@abstractmethod
def generate_tallies(self, materials, scores):
"""Use the C API to build tallies needed for reaction rates"""
@property
def nuclides(self):
"""List of nuclides with requested reaction rates"""
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
self._rate_tally.nuclides = nuclides
@abstractmethod
def get_material_rates(self, mat_id, nuc_index, react_index):
"""Return 2D array of [nuclide, reaction] reaction rates
Parameters
----------
mat_id : int
Unique ID for the requested material
nuc_index : list of str
Ordering of desired nuclides
react_index : list of str
Ordering of reactions
"""
def divide_by_adens(self, number):
"""Normalize reaction rates by number of nuclides
Acts on the current material examined by
:meth:`get_material_rates`
Parameters
----------
number : iterable of float
Number density [atoms/b-cm] of each nuclide tracked in the calculation.
Returns
-------
results : numpy.ndarray
Array of reactions rates of shape ``(n_nuclides, n_rxns)``
normalized by the number of nuclides
"""
mask = nonzero(number)
results = self._results_cache
for col in range(results.shape[1]):
results[mask, col] /= number[mask]
return results
class EnergyHelper(ABC):
"""Abstract class for obtaining energy produced
The ultimate goal of this helper is to provide instances of
:class:`openmc.deplete.Operator` with the total energy produced
in a transport simulation. This information, provided with the
power requested by the user and reaction rates from a
:class:`ReactionRateHelper` will scale reaction rates to the
correct values.
Attributes
----------
nuclides : list of str
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
energy : float
Total energy [J/s/source neutron] produced in a transport simulation.
Updated in the material iteration with :meth:`update`.
"""
def __init__(self):
self._nuclides = None
self._energy = 0.0
@property
def energy(self):
return self._energy * JOULE_PER_EV
def reset(self):
"""Reset energy produced prior to unpacking tallies"""
self._energy = 0.0
@abstractmethod
def prepare(self, chain_nucs, rate_index, materials):
"""Perform work needed to obtain energy produced
This method is called prior to the transport simulations
in :meth:`openmc.deplete.Operator.initial_condition`.
Parameters
----------
chain_nucs : list of str
All nuclides to be tracked in this problem
rate_index : dict of str to int
Mapping from nuclide name to index in the
`fission_rates` for :meth:`update`.
materials : list of str
All materials tracked on the operator helped by this
object. Should correspond to
:attr:`openmc.deplete.Operator.burnable_materials`
"""
def update(self, fission_rates, mat_index):
"""Update the energy produced
Parameters
----------
fission_rates : numpy.ndarray
fission reaction rate for each isotope in the specified
material. Should be ordered corresponding to initial
``rate_index`` used in :meth:`prepare`
mat_index : int
Index for the specific material in the list of all burnable
materials.
"""
@property
def nuclides(self):
"""List of nuclides with requested reaction rates"""
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides

164
openmc/deplete/chain.py
View file

@ -10,8 +10,10 @@ import math
import re
from collections import OrderedDict, defaultdict
from collections.abc import Mapping
from warnings import warn
from openmc.checkvalue import check_type
from openmc.checkvalue import check_type, check_less_than
from openmc.data import gnd_name, zam
# Try to use lxml if it is available. It preserves the order of attributes and
# provides a pretty-printer by default. If not available,
@ -451,3 +453,163 @@ class Chain(object):
matrix_dok = sp.dok_matrix((n, n))
dict.update(matrix_dok, matrix)
return matrix_dok.tocsr()
def get_capture_branches(self):
"""Return a dictionary with capture branching ratios
Returns
-------
capt :
nested dict of parent nuclide keys with capture targets and
branching ratios::
{"Am241": {"Am242": 0.91, "Am242_m1": 0.09}}
See Also
--------
:meth:`set_capture_branches`
"""
capt = {}
for nuclide in self.nuclides:
nuc_capt = {}
for rx in nuclide.reactions:
if rx.type == "(n,gamma)" and rx.branching_ratio != 1.0:
nuc_capt[rx.target] = rx.branching_ratio
if len(nuc_capt) > 0:
capt[nuclide.name] = nuc_capt
return capt
def set_capture_branches(self, branch_ratios, strict=True):
"""Set the capture branching ratios
To provide a buffer around floating point precisions,
the sum of all branching ratios from a single parent
cannot be greater than 1.00001.
Parameters
----------
branch_ratios : dict of {str: {str: float}}
Capture branching ratios to be inserted.
First layer keys are names of parent nuclides, e.g.
``"Am241"``. The capture branching ratios for these
parents will be modified. Corresponding values are
dictionaries of ``{target: branching_ratio}``
strict : bool
If this evalutes to ``True``, then all parents and
products must exist in the :class:`Chain`. A
:class:`KeyError` will be raised at the first
nuclide that does not exist. Otherwise, print
a warning message for missing parents and/or
products.
See Also
--------
:meth:`get_capture_branches`
"""
# Store some useful information through the validation stage
sums = {}
capt_ix_map = {}
grounds = {}
missing_parents = set()
missing_products = {}
no_capture = set()
# Check for validity before manipulation
for parent, sub in branch_ratios.items():
if parent not in self:
if strict:
raise KeyError(parent)
missing_parents.add(parent)
continue
# Make sure all products are present in the chain
prod_flag = False
for product in sub:
if product not in self:
if strict:
raise KeyError(product)
missing_products[parent] = product
prod_flag = True
break
if prod_flag:
continue
# Make sure this nuclide has capture reactions
indexes = []
for ix, rx in enumerate(self[parent].reactions):
if rx.type == "(n,gamma)":
indexes.append(ix)
if "_m" not in rx.target:
grounds[parent] = rx.target
if len(indexes) == 0:
if strict:
raise AttributeError(
"Nuclide {} does not have capture reactions in "
"this {}".format(parent, self.__class__.__name__))
no_capture.add(parent)
continue
capt_ix_map[parent] = indexes
this_sum = sum(sub.values())
check_less_than(parent + " ratios", this_sum, 1.00001)
sums[parent] = this_sum
if len(missing_parents) > 0:
warn("The following nuclides were not found in {}: {}".format(
self.__class__.__name__, ", ".join(sorted(missing_parents))))
if len(no_capture) > 0:
warn("The following nuclides did not have capture reactions: "
"{}".format(", ".join(sorted(no_capture))))
if len(missing_products) > 0:
tail = ("{} -> {}".format(k, v)
for k, v in sorted(missing_products.items()))
warn("The following products were not found in the {} and "
"parents were unmodified: \n{}".format(
self.__class__.__name__, ", ".join(tail)))
# Insert new ReactionTuples with updated branch ratios
for parent_name, capt_index in capt_ix_map.items():
parent = self[parent_name]
new_ratios = branch_ratios[parent_name]
capt_index = capt_ix_map[parent_name]
# Assume Q value is independent of target state
capt_Q = parent.reactions[capt_index[0]].Q
# Remove existing capture reactions
for ix in reversed(capt_index):
parent.reactions.pop(ix)
all_meta = True
for tgt, br in new_ratios.items():
all_meta = all_meta and ("_m" in tgt)
parent.reactions.append(ReactionTuple(
"(n,gamma)", tgt, capt_Q, br))
if all_meta and sums[parent_name] != 1.0:
ground_br = 1.0 - sums[parent_name]
ground_tgt = grounds.get(parent_name)
if ground_tgt is None:
pz, pa, pm = zam(parent_name)
ground_tgt = gnd_name(pz, pa + 1, 0)
new_ratios[ground_tgt] = ground_br
parent.reactions.append(ReactionTuple(
"(n,gamma)", ground_tgt, capt_Q, ground_br))

144
openmc/deplete/helpers.py Normal file
View file

@ -0,0 +1,144 @@
"""
Class for normalizing fission energy deposition
"""
from itertools import product
from numpy import dot, zeros
from openmc.capi import Tally, MaterialFilter
from .abc import ReactionRateHelper, EnergyHelper
# -------------------------------------
# Helpers for generating reaction rates
# -------------------------------------
class DirectReactionRateHelper(ReactionRateHelper):
"""Class that generates tallies for one-group rates
Parameters
----------
n_nucs : int
Number of burnable nuclides tracked by :class:`openmc.deplete.Operator`
n_react : int
Number of reactions tracked by :class:`openmc.deplete.Operator`
Attributes
----------
nuclides : list of str
All nuclides with desired reaction rates.
"""
def generate_tallies(self, materials, scores):
"""Produce one-group reaction rate tally
Uses the :mod:`openmc.capi` to generate a tally
of relevant reactions across all burnable materials.
Parameters
----------
materials : iterable of :class:`openmc.Material`
Burnable materials in the problem. Used to
construct a :class:`openmc.MaterialFilter`
scores : iterable of str
Reaction identifiers, e.g. ``"(n, fission)"``,
``"(n, gamma)"``, needed for the reaction rate tally.
"""
self._rate_tally = Tally()
self._rate_tally.scores = scores
self._rate_tally.filters = [MaterialFilter(materials)]
def get_material_rates(self, mat_id, nuc_index, react_index):
"""Return an array of reaction rates for a material
Parameters
----------
mat_id : int
Unique ID for the requested material
nuc_index : iterable of int
Index for each nuclide in :attr:`nuclides` in the
desired reaction rate matrix
react_index : iterable of int
Index for each reaction scored in the tally
Returns
-------
rates : numpy.ndarray
Array with shape ``(n_nuclides, n_rxns)`` with the
reaction rates in this material
"""
self._results_cache.fill(0.0)
full_tally_res = self._rate_tally.results[mat_id, :, 1]
for i_tally, (i_nuc, i_react) in enumerate(
product(nuc_index, react_index)):
self._results_cache[i_nuc, i_react] = full_tally_res[i_tally]
return self._results_cache
# ----------------------------
# Helpers for obtaining energy
# ----------------------------
class ChainFissionHelper(EnergyHelper):
"""Computes energy using fission Q values from depletion chain
Attributes
----------
nuclides : list of str
All nuclides with desired reaction rates. Ordered to be
consistent with :class:`openmc.deplete.Operator`
energy : float
Total energy [J/s/source neutron] produced in a transport simulation.
Updated in the material iteration with :meth:`update`.
"""
def __init__(self):
super().__init__()
self._fission_q_vector = None
def prepare(self, chain_nucs, rate_index, _materials):
"""Populate the fission Q value vector from a chain.
Parameters
----------
chain_nucs : iterable of :class:`openmc.deplete.Nuclide`
Nuclides used in this depletion chain. Do not need
to be ordered
rate_index : dict of str to int
Dictionary mapping names of nuclides, e.g. ``"U235"``,
to a corresponding index in the desired fission Q
vector.
_materials : list of str
Unused. Materials to be tracked for this helper.
"""
if (self._fission_q_vector is not None
and self._fission_q_vector.shape == (len(rate_index),)):
return
fission_qs = zeros(len(rate_index))
for nuclide in chain_nucs:
if nuclide.name in rate_index:
for rx in nuclide.reactions:
if rx.type == "fission":
fission_qs[rate_index[nuclide.name]] = rx.Q
break
self._fission_q_vector = fission_qs
def update(self, fission_rates, _mat_index):
"""Update energy produced with fission rates in a material
Parameters
----------
fission_rates : numpy.ndarray
fission reaction rate for each isotope in the specified
material. Should be ordered corresponding to initial
``rate_index`` used in :meth:`prepare`
_mat_index : int
index for the material requested. Unused, as identical
isotopes in all materials have the same Q value.
"""
self._energy += dot(fission_rates, self._fission_q_vector)

113
openmc/deplete/operator.py
View file

@ -20,11 +20,11 @@ from uncertainties import ufloat
import openmc
import openmc.capi
from openmc.data import JOULE_PER_EV
from . import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .reaction_rates import ReactionRates
from .helpers import DirectReactionRateHelper, ChainFissionHelper
def _distribute(items):
@ -77,6 +77,11 @@ class Operator(TransportOperator):
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``.
dilute_initial : float, optional
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
Defaults to 1.0e3.
Attributes
----------
@ -85,9 +90,9 @@ class Operator(TransportOperator):
settings : openmc.Settings
OpenMC settings object
dilute_initial : float
Initial atom density to add for nuclides that are zero in initial
condition to ensure they exist in the decay chain. Only done for
nuclides with reaction rates. Defaults to 1.0e3.
Initial atom density [atoms/cm^3] to add for nuclides that
are zero in initial condition to ensure they exist in the decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
@ -111,11 +116,11 @@ class Operator(TransportOperator):
Results from a previous depletion calculation
diff_burnable_mats : bool
Whether to differentiate burnable materials with multiple instances
"""
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False, fission_q=None):
super().__init__(chain_file, fission_q)
diff_burnable_mats=False, fission_q=None,
dilute_initial=1.0e3):
super().__init__(chain_file, fission_q, dilute_initial)
self.round_number = False
self.settings = settings
self.geometry = geometry
@ -139,6 +144,11 @@ class Operator(TransportOperator):
self.burnable_mats, volume, nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
# Generate map from local materials => material index
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
# Determine which nuclides have incident neutron data
self.nuclides_with_data = self._get_nuclides_with_data()
@ -153,6 +163,11 @@ class Operator(TransportOperator):
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
# Get classes to assist working with tallies
self._rate_helper = DirectReactionRateHelper(
self.reaction_rates.n_nuc, self.reaction_rates.n_react)
self._energy_helper = ChainFissionHelper()
def __call__(self, vec, power):
"""Runs a simulation.
@ -179,7 +194,8 @@ class Operator(TransportOperator):
# Update material compositions and tally nuclides
self._update_materials()
self._tally.nuclides = self._get_tally_nuclides()
self._rate_helper.nuclides = self._get_tally_nuclides()
self._energy_helper.nuclides = self._rate_helper.nuclides
# Run OpenMC
openmc.capi.reset()
@ -365,7 +381,11 @@ class Operator(TransportOperator):
openmc.capi.init(intracomm=comm)
# Generate tallies in memory
self._generate_tallies()
materials = [openmc.capi.materials[int(i)]
for i in self.burnable_mats]
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._energy_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc, materials)
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
@ -474,27 +494,6 @@ class Operator(TransportOperator):
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _generate_tallies(self):
"""Generates depletion tallies.
Using information from the depletion chain as well as the nuclides
currently in the problem, this function automatically generates a
tally.xml for the simulation.
"""
# Create tallies for depleting regions
materials = [openmc.capi.materials[int(i)]
for i in self.burnable_mats]
mat_filter = openmc.capi.MaterialFilter(materials)
# Set up a tally that has a material filter covering each depletable
# material and scores corresponding to all reactions that cause
# transmutation. The nuclides for the tally are set later when eval() is
# called.
self._tally = openmc.capi.Tally()
self._tally.scores = self.chain.reactions
self._tally.filters = [mat_filter]
def _unpack_tallies_and_normalize(self, power):
"""Unpack tallies from OpenMC and return an operator result
@ -515,78 +514,54 @@ class Operator(TransportOperator):
"""
rates = self.reaction_rates
rates[:, :, :] = 0.0
rates.fill(0.0)
# Get k and uncertainty
k_combined = ufloat(*openmc.capi.keff())
# Extract tally bins
materials = self.burnable_mats
nuclides = self._tally.nuclides
nuclides = self._rate_helper.nuclides
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Compute fission power
# TODO : improve this calculation
# Keep track of energy produced from all reactions in eV per source
# particle
energy = 0.0
self._energy_helper.reset()
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers
fission_Q = np.zeros(rates.n_nuc)
rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
number = np.zeros(rates.n_nuc)
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx["fission"]
for nuclide in self.chain.nuclides:
if nuclide.name in rates.index_nuc:
for rx in nuclide.reactions:
if rx.type == 'fission':
ind = rates.index_nuc[nuclide.name]
fission_Q[ind] = rx.Q
break
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
slab = materials.index(mat)
# Get material results hyperslab
results = self._tally.results[slab, :, 1]
mat_index = self._mat_index_map[mat]
# Zero out reaction rates and nuclide numbers
rates_expanded[:] = 0.0
number[:] = 0.0
number.fill(0.0)
# Expand into our memory layout
j = 0
# Get new number densities
for nuc, i_nuc_results in zip(nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
for react in react_ind:
rates_expanded[i_nuc_results, react] = results[j]
j += 1
tally_rates = self._rate_helper.get_material_rates(
mat_index, nuc_ind, react_ind)
# Accumulate energy from fission
energy += np.dot(rates_expanded[:, fission_ind], fission_Q)
self._energy_helper.update(tally_rates[:, fission_ind], mat_index)
# Divide by total number and store
for i_nuc_results in nuc_ind:
if number[i_nuc_results] != 0.0:
for react in react_ind:
rates_expanded[i_nuc_results, react] /= number[i_nuc_results]
rates[i, :, :] = rates_expanded
rates[i] = self._rate_helper.divide_by_adens(number)
# Reduce energy produced from all processes
energy = comm.allreduce(energy)
# Determine power in eV/s
power /= JOULE_PER_EV
# J / s / source neutron
energy = comm.allreduce(self._energy_helper.energy)
# Scale reaction rates to obtain units of reactions/sec
rates *= power / energy

18
openmc/deplete/results_list.py
View file

@ -26,7 +26,15 @@ class ResultsList(list):
self.append(Results.from_hdf5(fh, i))
def get_atoms(self, mat, nuc):
"""Get nuclide concentration over time from a single material
"""Get number of nuclides over time from a single material
.. note::
Initial values for some isotopes that do not appear in
initial concentrations may be non-zero, depending on the
value of :class:`openmc.deplete.Operator` ``dilute_initial``.
The :class:`openmc.deplete.Operator` adds isotopes according
to this setting, which can be set to zero.
Parameters
----------
@ -56,6 +64,14 @@ class ResultsList(list):
def get_reaction_rate(self, mat, nuc, rx):
"""Get reaction rate in a single material/nuclide over time
.. note::
Initial values for some isotopes that do not appear in
initial concentrations may be non-zero, depending on the
value of :class:`openmc.deplete.Operator` ``dilute_initial``
The :class:`openmc.deplete.Operator` adds isotopes according
to this setting, which can be set to zero.
Parameters
----------
mat : str

30
openmc/filter.py
View file

@ -32,7 +32,8 @@ _CURRENT_NAMES = (
'z-min out', 'z-min in', 'z-max out', 'z-max in'
)
_PARTICLE_IDS = {'neutron': 1, 'photon': 2, 'electron': 3, 'positron': 4}
_PARTICLES = {'neutron', 'photon', 'electron', 'positron'}
class FilterMeta(ABCMeta):
def __new__(cls, name, bases, namespace, **kwargs):
@ -547,10 +548,9 @@ class ParticleFilter(Filter):
Parameters
----------
bins : str, int, or iterable of Integral
The Particles to tally. Either str with particle type or their
ID numbers can be used ('neutron' = 1, 'photon' = 2, 'electron' = 3,
'positron' = 4).
bins : str, or iterable of str
The particles to tally represented as strings ('neutron', 'photon',
'electron', 'positron').
filter_id : int
Unique identifier for the filter
@ -571,16 +571,22 @@ class ParticleFilter(Filter):
@bins.setter
def bins(self, bins):
bins = np.atleast_1d(bins)
cv.check_iterable_type('filter bins', bins, (Integral, str))
cv.check_iterable_type('filter bins', bins, str)
for edge in bins:
if isinstance(edge, Integral):
cv.check_value('filter bin', edge, _PARTICLE_IDS.values())
else:
cv.check_value('filter bin', edge, _PARTICLE_IDS.keys())
bins = np.atleast_1d([b if isinstance(b, Integral) else _PARTICLE_IDS[b]
for b in bins])
cv.check_value('filter bin', edge, _PARTICLES)
self._bins = bins
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'][()].decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'][()].decode() + " instead")
particles = [b.decode() for b in group['bins'][()]]
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
return cls(particles, filter_id=filter_id)
class MeshFilter(Filter):
"""Bins tally event locations onto a regular, rectangular mesh.

11
openmc/material.py
View file

@ -287,6 +287,8 @@ class Material(IDManagerMixin):
material.depletable = bool(group.attrs['depletable'])
if 'volume' in group.attrs:
material.volume = group.attrs['volume']
if "temperature" in group.attrs:
material.temperature = group.attrs["temperature"]
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in group:
@ -841,8 +843,7 @@ class Material(IDManagerMixin):
# Create temperature XML subelement
if self.temperature is not None:
subelement = ET.SubElement(element, "temperature")
subelement.text = str(self.temperature)
element.set("temperature", str(self.temperature))
# Create density XML subelement
if self._density is not None or self._density_units == 'sum':
@ -931,8 +932,10 @@ class Material(IDManagerMixin):
mat_id = int(elem.get('id'))
mat = cls(mat_id)
mat.name = elem.get('name')
if 'temperature' in elem.attrib:
mat.temperature = float(elem.get('temperature'))
if "temperature" in elem.attrib:
mat.temperature = float(elem.get("temperature"))
if 'volume' in elem.attrib:
mat.volume = float(elem.get('volume'))
mat.depletable = bool(elem.get('depletable'))

2
openmc/mgxs/mgxs.py
View file

@ -2751,7 +2751,7 @@ class TransportXS(MGXS):
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
filter_bins=[('P1',)],
squeeze=True)
p1_tally.scores = ['scatter-1']
p1_tally._scores = ['scatter-1']
self._rxn_rate_tally = self.tallies['total'] - p1_tally
self._rxn_rate_tally.sparse = self.sparse

156
openmc/model/funcs.py
View file

@ -1,12 +1,16 @@
from collections import OrderedDict
from collections.abc import Iterable
from math import sqrt
from numbers import Real
from functools import partial
from warnings import warn
from operator import attrgetter
from openmc import XPlane, YPlane, Plane, ZCylinder, Quadric
from openmc.checkvalue import check_type, check_value
from openmc import (
XPlane, YPlane, Plane, ZCylinder, Quadric, Cylinder, XCylinder,
YCylinder, Material, Universe, Cell)
from openmc.checkvalue import (
check_type, check_value, check_length, check_less_than,
check_iterable_type)
import openmc.data
@ -395,9 +399,9 @@ def cylinder_from_points(p1, p2, r, **kwargs):
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
cx = y1*z2 + y2*z1
cy = -(x1*z2 + x2*z1)
cz = x1*y2 + x2*y1
cx = y1*z2 - y2*z1
cy = x2*z1 - x1*z2
cz = x1*y2 - x2*y1
# Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the
# cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|.
@ -409,12 +413,12 @@ def cylinder_from_points(p1, p2, r, **kwargs):
kwargs['d'] = -2*dx*dy
kwargs['e'] = -2*dy*dz
kwargs['f'] = -2*dx*dz
kwargs['g'] = cy*dz - cz*dy
kwargs['h'] = cz*dx - cx*dz
kwargs['j'] = cx*dy - cy*dx
kwargs['k'] = -(dx*dx + dy*dy + dz*dz)*r*r
kwargs['g'] = 2*(cy*dz - cz*dy)
kwargs['h'] = 2*(cz*dx - cx*dz)
kwargs['j'] = 2*(cx*dy - cy*dx)
kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r
return openmc.Quadric(**kwargs)
return Quadric(**kwargs)
def subdivide(surfaces):
@ -442,3 +446,133 @@ def subdivide(surfaces):
regions.append(+s0 & -s1)
regions.append(+surfaces[-1])
return regions
def pin(surfaces, items, subdivisions=None, divide_vols=True,
**kwargs):
"""Convenience function for building a fuel pin
Parameters
----------
surfaces : iterable of :class:`openmc.Cylinder`
Cylinders used to define boundaries
between items. All cylinders must be
concentric and of the same orientation, e.g.
all :class:`openmc.ZCylinder`
items : iterable
Objects to go between ``surfaces``. These can be anything
that can fill a :class:`openmc.Cell`, including
:class:`openmc.Material`, or other :class:`openmc.Universe`
objects. There must be one more item than surfaces,
which will span all space outside the final ring.
subdivisions : None or dict of int to int
Dictionary describing which rings to subdivide and how
many times. Keys are indexes of the annular rings
to be divided. Will construct equal area rings
divide_vols : bool
If this evaluates to ``True``, then volumes of subdivided
:class:`openmc.Material` instances will also be divided by the
number of divisions. Otherwise the volume of the
original material will not be modified before subdivision
kwargs:
Additional key-word arguments to be passed to
:class:`openmc.Universe`, like ``name="Fuel pin"``
Returns
-------
:class:`openmc.Universe`
Universe of concentric cylinders filled with the desired
items
"""
if "cells" in kwargs:
raise SyntaxError(
"Cells will be set by this function, not from input arguments.")
check_type("items", items, Iterable)
check_length("surfaces", surfaces, len(items) - 1, len(items) - 1)
# Check that all surfaces are of similar orientation
check_type("surface", surfaces[0], Cylinder)
surf_type = type(surfaces[0])
check_iterable_type("surfaces", surfaces[1:], surf_type)
# Check for increasing radii and equal centers
if surf_type is ZCylinder:
center_getter = attrgetter("x0", "y0")
elif surf_type is YCylinder:
center_getter = attrgetter("x0", "z0")
elif surf_type is XCylinder:
center_getter = attrgetter("z0", "y0")
else:
raise TypeError(
"Not configured to interpret {} surfaces".format(
surf_type.__name__))
centers = set()
prev_rad = 0
for ix, surf in enumerate(surfaces):
cur_rad = surf.r
if cur_rad <= prev_rad:
raise ValueError(
"Surfaces do not appear to be increasing in radius. "
"Surface {} at index {} has radius {:7.3e} compared to "
"previous radius of {:7.5e}".format(
surf.id, ix, cur_rad, prev_rad))
prev_rad = cur_rad
centers.add(center_getter(surf))
if len(centers) > 1:
raise ValueError(
"Surfaces do not appear to be concentric. The following "
"centers were found: {}".format(centers))
if subdivisions is not None:
check_length("subdivisions", subdivisions, 1, len(surfaces))
orig_indexes = list(subdivisions.keys())
check_iterable_type("ring indexes", orig_indexes, int)
check_iterable_type(
"number of divisions", list(subdivisions.values()), int)
for ix in orig_indexes:
if ix < 0:
subdivisions[len(surfaces) + ix] = subdivisions.pop(ix)
# Dissallow subdivision on outer most, infinite region
check_less_than(
"outer ring", max(subdivisions), len(surfaces), equality=True)
# ensure ability to concatenate
if not isinstance(items, list):
items = list(items)
if not isinstance(surfaces, list):
surfaces = list(surfaces)
# generate equal area divisions
# Adding N - 1 new regions
# N - 2 surfaces are made
# Original cell is not removed, but not occupies last ring
for ring_index in reversed(sorted(subdivisions.keys())):
nr = subdivisions[ring_index]
new_surfs = []
lower_rad = 0.0 if ring_index == 0 else surfaces[ring_index - 1].r
upper_rad = surfaces[ring_index].r
area_term = (upper_rad ** 2 - lower_rad ** 2) / nr
for new_index in range(nr - 1):
lower_rad = sqrt(area_term + lower_rad ** 2)
new_surfs.append(surf_type(r=lower_rad))
surfaces = (
surfaces[:ring_index] + new_surfs + surfaces[ring_index:])
filler = items[ring_index]
if (divide_vols and hasattr(filler, "volume")
and filler.volume is not None):
filler.volume /= nr
items[ring_index:ring_index] = [
filler.clone() for _i in range(nr - 1)]
# Build the universe
regions = subdivide(surfaces)
cells = [Cell(fill=f, region=r) for r, f in zip(regions, items)]
return Universe(cells=cells, **kwargs)

242
src/cell.cpp
View file

@ -1,11 +1,12 @@
#include "openmc/cell.h"
#include <cctype>
#include <cmath>
#include <sstream>
#include <set>
#include <string>
#include <cctype>
#include <gsl/gsl>
#include "openmc/capi.h"
#include "openmc/constants.h"
@ -118,7 +119,6 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
if (token < OP_UNION) {
// If token is not an operator, add it to output
rpn.push_back(token);
} else if (token < OP_RIGHT_PAREN) {
// Regular operators union, intersection, complement
while (stack.size() > 0) {
@ -158,7 +158,6 @@ generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
<< cell_id;
fatal_error(err_msg);
}
rpn.push_back(stack.back());
stack.pop_back();
}
@ -208,10 +207,56 @@ Universe::to_hdf5(hid_t universes_group) const
close_group(group);
}
BoundingBox Universe::bounding_box() const {
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
return {};
} else {
for (const auto& cell : cells_) {
auto& c = model::cells[cell];
bbox |= c->bounding_box();
}
}
return bbox;
}
//==============================================================================
// Cell implementation
//==============================================================================
double
Cell::temperature(int32_t instance) const
{
if (sqrtkT_.size() < 1) {
throw std::runtime_error{"Cell temperature has not yet been set."};
}
if (instance >= 0) {
double sqrtkT = sqrtkT_.size() == 1 ?
sqrtkT_.at(0) :
sqrtkT_.at(instance);
return sqrtkT * sqrtkT / K_BOLTZMANN;
} else {
return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN;
}
}
void
Cell::set_temperature(double T, int32_t instance)
{
if (instance >= 0) {
sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
} else {
for (auto& T_ : sqrtkT_) {
T_ = std::sqrt(K_BOLTZMANN * T);
}
}
}
//==============================================================================
// CSGCell implementation
//==============================================================================
CSGCell::CSGCell() {} // empty constructor
CSGCell::CSGCell(pugi::xml_node cell_node)
@ -534,6 +579,95 @@ CSGCell::to_hdf5(hid_t cell_group) const
close_group(group);
}
BoundingBox CSGCell::bounding_box_simple() const {
BoundingBox bbox;
for (int32_t token : rpn_) {
bbox &= model::surfaces[abs(token)-1]->bounding_box(token > 0);
}
return bbox;
}
void CSGCell::apply_demorgan(std::vector<int32_t>& rpn) {
for (auto& token : rpn) {
if (token < OP_UNION) { token *= -1; }
else if (token == OP_UNION) { token = OP_INTERSECTION; }
else if (token == OP_INTERSECTION) { token = OP_UNION; }
}
}
BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
// if the last operator is a complement op, there is no
// sub-region that the complement connects to. This indicates
// that the entire region is a complement and we can apply
// De Morgan's laws immediately
if (rpn.back() == OP_COMPLEMENT) {
rpn.pop_back();
apply_demorgan(rpn);
}
// reverse the rpn to make popping easier
std::reverse(rpn.begin(), rpn.end());
BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0);
rpn.pop_back();
while (rpn.size()) {
// move through the rpn in twos
int32_t one = rpn.back(); rpn.pop_back();
int32_t two = rpn.back(); rpn.pop_back();
// the first token should always be a surface
Expects(one < OP_UNION);
if (two >= OP_UNION) {
if (two == OP_UNION) {
current |= model::surfaces[abs(one)-1]->bounding_box(one > 0);
} else if (two == OP_INTERSECTION) {
current &= model::surfaces[abs(one)-1]->bounding_box(one > 0);
}
} else {
// two surfaces in a row (left parenthesis),
// create sub-rpn for region in parenthesis
std::vector<int32_t> subrpn;
subrpn.push_back(one);
subrpn.push_back(two);
// add until last two tokens in the sub-rpn are operators
// (indicates a right parenthesis)
while (!((subrpn.back() >= OP_UNION) && (*(subrpn.rbegin() + 1) >= OP_UNION))) {
subrpn.push_back(rpn.back());
rpn.pop_back();
}
// handle complement case using De Morgan's laws
if (subrpn.back() == OP_COMPLEMENT) {
subrpn.pop_back();
apply_demorgan(subrpn);
subrpn.push_back(rpn.back());
rpn.pop_back();
}
// save the last operator, tells us how to combine this region
// with our current bounding box
int32_t op = subrpn.back(); subrpn.pop_back();
// get bounding box for the subrpn
BoundingBox sub_box = bounding_box_complex(subrpn);
// combine the sub-rpn bounding box with our current cell box
if (op == OP_UNION) {
current |= sub_box;
} else if (op == OP_INTERSECTION) {
current &= sub_box;
}
}
}
return current;
}
BoundingBox CSGCell::bounding_box() const {
return simple_ ? bounding_box_simple() : bounding_box_complex(rpn_);
}
//==============================================================================
bool
@ -658,6 +792,16 @@ bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
void DAGCell::to_hdf5(hid_t group_id) const { return; }
BoundingBox DAGCell::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(vol, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
#endif
//==============================================================================
@ -917,27 +1061,18 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
extern "C" int
openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
{
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
if (instance) {
if (*instance >= 0 && *instance < c.sqrtkT_.size()) {
c.sqrtkT_[*instance] = std::sqrt(K_BOLTZMANN * T);
} else {
strcpy(openmc_err_msg, "Distribcell instance is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
} else {
for (auto& T_ : c.sqrtkT_) {
T_ = std::sqrt(K_BOLTZMANN * T);
}
}
} else {
if (index < 0 || index >= model::cells.size()) {
strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
int32_t instance_index = instance ? *instance : -1;
try {
model::cells[index]->set_temperature(T, instance_index);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
return 0;
}
@ -949,28 +1084,65 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
return OPENMC_E_OUT_OF_BOUNDS;
}
Cell& c {*model::cells[index]};
if (c.sqrtkT_.size() < 1) {
strcpy(openmc_err_msg, "Cell temperature has not yet been set.");
int32_t instance_index = instance ? *instance : -1;
try {
*T = model::cells[index]->temperature(instance_index);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
return 0;
}
if (instance) {
if (*instance >= 0 && *instance < c.n_instances_) {
double sqrtkT = c.sqrtkT_.size() == 1 ? c.sqrtkT_[0] : c.sqrtkT_[*instance];
*T = sqrtkT * sqrtkT / K_BOLTZMANN;
} else {
strcpy(openmc_err_msg, "Distribcell instance is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
} else {
*T = c.sqrtkT_[0] * c.sqrtkT_[0] / K_BOLTZMANN;
}
//! Get the bounding box of a cell
extern "C" int
openmc_cell_bounding_box(const int32_t index, double* llc, double* urc) {
BoundingBox bbox;
const auto& c = model::cells[index];
bbox = c->bounding_box();
// set lower left corner values
llc[0] = bbox.xmin;
llc[1] = bbox.ymin;
llc[2] = bbox.zmin;
// set upper right corner values
urc[0] = bbox.xmax;
urc[1] = bbox.ymax;
urc[2] = bbox.zmax;
return 0;
}
//! Get the name of a cell
extern "C" int
openmc_cell_get_name(int32_t index, const char** name) {
if (index < 0 || index >= model::cells.size()) {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*name = model::cells[index]->name().data();
return 0;
}
//! Set the name of a cell
extern "C" int
openmc_cell_set_name(int32_t index, const char* name) {
if (index < 0 || index >= model::cells.size()) {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::cells[index]->set_name(name);
return 0;
}
//! Return the index in the cells array of a cell with a given ID
extern "C" int
openmc_get_cell_index(int32_t id, int32_t* index)

60
src/dagmc.cpp
View file

@ -35,10 +35,9 @@ const bool dagmc_enabled = false;
#ifdef DAGMC
const std::string DAGMC_FILENAME = "dagmc.h5m";
namespace openmc {
const std::string DAGMC_FILENAME = "dagmc.h5m";
namespace simulation {
@ -54,8 +53,16 @@ moab::DagMC* DAG;
} // namespace model
void check_dagmc_file() {
std::string filename = settings::path_input + DAGMC_FILENAME;
if (!file_exists(filename)) {
fatal_error("Geometry DAGMC file '" + filename + "' does not exist!");
}
}
bool get_uwuw_materials_xml(std::string& s) {
UWUW uwuw(DAGMC_FILENAME.c_str());
check_dagmc_file();
UWUW uwuw((settings::path_input + DAGMC_FILENAME).c_str());
std::stringstream ss;
bool uwuw_mats_present = false;
@ -133,7 +140,7 @@ void legacy_assign_material(const std::string& mat_string, DAGCell* c)
}
if (settings::verbosity >= 10) {
Material* m = model::materials[model::material_map[c->material_[0]]].get();
const auto& m = model::materials[model::material_map.at(c->material_[0])];
std::stringstream msg;
msg << "DAGMC material " << mat_string << " was assigned";
if (mat_found_by_name) {
@ -147,14 +154,18 @@ void legacy_assign_material(const std::string& mat_string, DAGCell* c)
void load_dagmc_geometry()
{
check_dagmc_file();
if (!model::DAG) {
model::DAG = new moab::DagMC();
}
std::string filename = settings::path_input + DAGMC_FILENAME;
// --- Materials ---
// create uwuw instance
UWUW uwuw(DAGMC_FILENAME.c_str());
UWUW uwuw(filename.c_str());
// check for uwuw material definitions
bool using_uwuw = !uwuw.material_library.empty();
@ -167,7 +178,7 @@ void load_dagmc_geometry()
int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC runs
// load the DAGMC geometry
moab::ErrorCode rval = model::DAG->load_file(DAGMC_FILENAME.c_str());
moab::ErrorCode rval = model::DAG->load_file(filename.c_str());
MB_CHK_ERR_CONT(rval);
// initialize acceleration data structures
@ -215,17 +226,6 @@ void load_dagmc_geometry()
model::universes[it->second]->cells_.push_back(i);
}
// check for temperature assignment
std::string temp_value;
if (model::DAG->has_prop(vol_handle, "temp")) {
rval = model::DAG->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default));
}
// MATERIALS
if (model::DAG->is_implicit_complement(vol_handle)) {
@ -295,6 +295,26 @@ void load_dagmc_geometry()
legacy_assign_material(mat_value, c);
}
}
// check for temperature assignment
std::string temp_value;
// no temperature if void
if (c->material_[0] == MATERIAL_VOID) continue;
// assign cell temperature
const auto& mat = model::materials[model::material_map.at(c->material_[0])];
if (model::DAG->has_prop(vol_handle, "temp")) {
rval = model::DAG->prop_value(vol_handle, "temp", temp_value);
MB_CHK_ERR_CONT(rval);
double temp = std::stod(temp_value);
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * temp));
} else if (mat->temperature_ > 0.0) {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature_));
} else {
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * settings::temperature_default));
}
}
// allocate the cell overlap count if necessary
@ -369,11 +389,7 @@ void load_dagmc_geometry()
void read_geometry_dagmc()
{
// Check if dagmc.h5m exists
std::string filename = settings::path_input + "dagmc.h5m";
if (!file_exists(filename)) {
fatal_error("Geometry DAGMC file '" + filename + "' does not exist!");
}
check_dagmc_file();
write_message("Reading DAGMC geometry...", 5);
load_dagmc_geometry();

18
src/geometry.cpp
View file

@ -9,6 +9,7 @@
#include "openmc/lattice.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/surface.h"
@ -475,4 +476,21 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
return 0;
}
extern "C" int openmc_global_bounding_box(double* llc, double* urc) {
auto bbox = model::universes.at(model::root_universe)->bounding_box();
// set lower left corner values
llc[0] = bbox.xmin;
llc[1] = bbox.ymin;
llc[2] = bbox.zmin;
// set upper right corner values
urc[0] = bbox.xmax;
urc[1] = bbox.ymax;
urc[2] = bbox.zmax;
return 0;
}
} // namespace openmc

2
src/geometry_aux.cpp
View file

@ -309,7 +309,7 @@ prepare_distribcell()
for (auto& filt : model::tally_filters) {
auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
if (distrib_filt) {
distribcells.insert(distrib_filt->cell_);
distribcells.insert(distrib_filt->cell());
}
}

282
src/material.cpp
View file

@ -47,9 +47,10 @@ std::unordered_map<int32_t, int32_t> material_map;
//==============================================================================
Material::Material(pugi::xml_node node)
: index_{model::materials.size()}
{
if (check_for_node(node, "id")) {
id_ = std::stoi(get_node_value(node, "id"));
this->set_id(std::stoi(get_node_value(node, "id")));
} else {
fatal_error("Must specify id of material in materials XML file.");
}
@ -331,6 +332,11 @@ Material::Material(pugi::xml_node node)
}
}
Material::~Material()
{
model::material_map.erase(id_);
}
void Material::finalize()
{
// Set fissionable if any nuclide is fissionable
@ -849,11 +855,41 @@ void Material::calculate_photon_xs(Particle& p) const
}
}
int Material::set_density(double density, std::string units)
void Material::set_id(int32_t id)
{
Expects(id >= -1);
// Clear entry in material map if an ID was already assigned before
if (id_ != -1) {
model::material_map.erase(id_);
id_ = -1;
}
// Make sure no other material has same ID
if (model::material_map.find(id) != model::material_map.end()) {
throw std::runtime_error{"Two materials have the same ID: " + std::to_string(id)};
}
// If no ID specified, auto-assign next ID in sequence
if (id == -1) {
id = 0;
for (const auto& f : model::materials) {
id = std::max(id, f->id_);
}
++id;
}
// Update ID and entry in material map
id_ = id;
model::material_map[id] = index_;
}
void Material::set_density(double density, gsl::cstring_span units)
{
Expects(density >= 0.0);
if (nuclide_.empty()) {
set_errmsg("No nuclides exist in material yet.");
return OPENMC_E_ALLOCATE;
throw std::runtime_error{"No nuclides exist in material yet."};
}
if (units == "atom/b-cm") {
@ -884,10 +920,51 @@ int Material::set_density(double density, std::string units)
density_ *= f;
atom_density_ *= f;
} else {
set_errmsg("Invalid units '" + units + "' specified.");
return OPENMC_E_INVALID_ARGUMENT;
throw std::invalid_argument{"Invalid units '" + std::string(units.data())
+ "' specified."};
}
return 0;
}
void Material::set_densities(const std::vector<std::string>& name,
const std::vector<double>& density)
{
auto n = name.size();
Expects(n > 0);
Expects(n == density.size());
if (n != nuclide_.size()) {
nuclide_.resize(n);
atom_density_ = xt::zeros<double>({n});
}
double sum_density = 0.0;
for (gsl::index i = 0; i < n; ++i) {
const auto& nuc {name[i]};
if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) {
int err = openmc_load_nuclide(nuc.c_str());
if (err < 0) throw std::runtime_error{openmc_err_msg};
}
nuclide_[i] = data::nuclide_map.at(nuc);
Expects(density[i] > 0.0);
atom_density_(i) = density[i];
sum_density += density[i];
}
// Set total density to the sum of the vector
this->set_density(sum_density, "atom/b-cm");
// Assign S(a,b) tables
this->init_thermal();
}
double Material::volume() const
{
if (volume_ < 0.0) {
throw std::runtime_error{"Volume for material with ID="
+ std::to_string(id_) + " not set."};
}
return volume_;
}
void Material::to_hdf5(hid_t group) const
@ -898,6 +975,9 @@ void Material::to_hdf5(hid_t group) const
if (volume_ > 0.0) {
write_attribute(material_group, "volume", volume_);
}
if (temperature_ > 0.0) {
write_attribute(material_group, "temperature", temperature_);
}
write_dataset(material_group, "name", name_);
write_dataset(material_group, "atom_density", density_);
@ -945,6 +1025,42 @@ void Material::to_hdf5(hid_t group) const
close_group(material_group);
}
void Material::add_nuclide(const std::string& name, double density)
{
// Check if nuclide is already in material
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
if (data::nuclides[i_nuc]->name_ == name) {
double awr = data::nuclides[i_nuc]->awr_;
density_ += density - atom_density_(i);
density_gpcc_ += (density - atom_density_(i))
* awr * MASS_NEUTRON / N_AVOGADRO;
atom_density_(i) = density;
return;
}
}
// If nuclide wasn't found, extend nuclide/density arrays
int err = openmc_load_nuclide(name.c_str());
if (err < 0) throw std::runtime_error{openmc_err_msg};
// Append new nuclide/density
int i_nuc = data::nuclide_map[name];
nuclide_.push_back(i_nuc);
auto n = nuclide_.size();
// Create copy of atom_density_ array with one extra entry
xt::xtensor<double, 1> atom_density = xt::zeros<double>({n});
xt::view(atom_density, xt::range(0, n-1)) = atom_density_;
atom_density(n-1) = density;
atom_density_ = atom_density;
density_ += density;
density_gpcc_ += density * data::nuclides[i_nuc]->awr_
* MASS_NEUTRON / N_AVOGADRO;
}
//==============================================================================
// Non-method functions
//==============================================================================
@ -1098,19 +1214,6 @@ void read_materials_xml()
model::materials.push_back(std::make_unique<Material>(material_node));
}
model::materials.shrink_to_fit();
// Populate the material map.
for (int i = 0; i < model::materials.size(); i++) {
int32_t mid = model::materials[i]->id_;
auto search = model::material_map.find(mid);
if (search == model::material_map.end()) {
model::material_map[mid] = i;
} else {
std::stringstream err_msg;
err_msg << "Two or more materials use the same unique ID: " << mid;
fatal_error(err_msg);
}
}
}
void free_memory_material()
@ -1141,40 +1244,10 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density)
{
int err = 0;
if (index >= 0 && index < model::materials.size()) {
auto& m = model::materials[index];
// Check if nuclide is already in material
for (int i = 0; i < m->nuclide_.size(); ++i) {
int i_nuc = m->nuclide_[i];
if (data::nuclides[i_nuc]->name_ == name) {
double awr = data::nuclides[i_nuc]->awr_;
m->density_ += density - m->atom_density_(i);
m->density_gpcc_ += (density - m->atom_density_(i))
* awr * MASS_NEUTRON / N_AVOGADRO;
m->atom_density_(i) = density;
return 0;
}
}
// If nuclide wasn't found, extend nuclide/density arrays
err = openmc_load_nuclide(name);
if (err == 0) {
// Append new nuclide/density
int i_nuc = data::nuclide_map[name];
m->nuclide_.push_back(i_nuc);
auto n = m->nuclide_.size();
// Create copy of atom_density_ array with one extra entry
xt::xtensor<double, 1> atom_density = xt::zeros<double>({n});
xt::view(atom_density, xt::range(0, n-1)) = m->atom_density_;
atom_density(n-1) = density;
m->atom_density_ = atom_density;
m->density_ += density;
m->density_gpcc_ += density * data::nuclides[i_nuc]->awr_
* MASS_NEUTRON / N_AVOGADRO;
try {
model::materials[index]->add_nuclide(name, density);
} catch (const std::runtime_error& e) {
return OPENMC_E_DATA;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
@ -1184,14 +1257,14 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density)
}
extern "C" int
openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n)
openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n)
{
if (index >= 0 && index < model::materials.size()) {
auto& mat = model::materials[index];
if (!mat->nuclide_.empty()) {
*nuclides = mat->nuclide_.data();
*densities = mat->atom_density_.data();
*n = mat->nuclide_.size();
if (!mat->nuclides().empty()) {
*nuclides = mat->nuclides().data();
*densities = mat->densities().data();
*n = mat->nuclides().size();
return 0;
} else {
set_errmsg("Material atom density array has not been allocated.");
@ -1208,7 +1281,7 @@ openmc_material_get_density(int32_t index, double* density)
{
if (index >= 0 && index < model::materials.size()) {
auto& mat = model::materials[index];
*density = mat->density_gpcc_;
*density = mat->density_gpcc();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
@ -1220,7 +1293,7 @@ extern "C" int
openmc_material_get_fissionable(int32_t index, bool* fissionable)
{
if (index >= 0 && index < model::materials.size()) {
*fissionable = model::materials[index]->fissionable_;
*fissionable = model::materials[index]->fissionable();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
@ -1232,7 +1305,7 @@ extern "C" int
openmc_material_get_id(int32_t index, int32_t* id)
{
if (index >= 0 && index < model::materials.size()) {
*id = model::materials[index]->id_;
*id = model::materials[index]->id();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
@ -1244,16 +1317,13 @@ extern "C" int
openmc_material_get_volume(int32_t index, double* volume)
{
if (index >= 0 && index < model::materials.size()) {
auto& m = model::materials[index];
if (m->volume_ >= 0.0) {
*volume = m->volume_;
return 0;
} else {
std::stringstream msg;
msg << "Volume for material with ID=" << m->id_ << " not set.";
set_errmsg(msg);
try {
*volume = model::materials[index]->volume();
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
@ -1264,59 +1334,75 @@ extern "C" int
openmc_material_set_density(int32_t index, double density, const char* units)
{
if (index >= 0 && index < model::materials.size()) {
return model::materials[index]->set_density(density, units);
try {
model::materials[index]->set_density(density, units);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int
openmc_material_set_densities(int32_t index, int n, const char** name, const double* density)
{
if (index >= 0 && index < model::materials.size()) {
auto& mat {model::materials[index]};
if (n != mat->nuclide_.size()) {
mat->nuclide_.resize(n);
mat->atom_density_ = xt::zeros<double>({n});
try {
model::materials[index]->set_densities({name, name + n}, {density, density + n});
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
double sum_density = 0.0;
for (int i = 0; i < n; ++i) {
std::string nuc {name[i]};
if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) {
int err = openmc_load_nuclide(nuc.c_str());
if (err < 0) return err;
}
mat->nuclide_[i] = data::nuclide_map[nuc];
mat->atom_density_(i) = density[i];
sum_density += density[i];
}
// Set total density to the sum of the vector
int err = mat->set_density(sum_density, "atom/b-cm");
// Assign S(a,b) tables
mat->init_thermal();
return err;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int
openmc_material_set_id(int32_t index, int32_t id)
{
if (index >= 0 && index < model::materials.size()) {
model::materials[index]->id_ = id;
model::material_map[id] = index;
return 0;
try {
model::materials.at(index)->set_id(id);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int
openmc_material_get_name(int32_t index, const char** name) {
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*name = model::materials[index]->name().data();
return 0;
}
extern "C" int
openmc_material_set_name(int32_t index, const char* name) {
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::materials[index]->set_name(name);
return 0;
}
extern "C" int

2
src/mgxs.cpp
View file

@ -522,7 +522,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs_t->decay_rate(a, *dg + 1);
val = xs_t->decay_rate(a, *dg);
} else {
val = xs_t->decay_rate(a, 0);
}

20
src/mgxs_interface.cpp
View file

@ -263,21 +263,13 @@ get_nuclide_xs(int index, int xstype, int gin, const int* gout,
{
int gout_c;
const int* gout_c_p;
int dg_c;
const int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return data::nuclides_MG[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
return data::nuclides_MG[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
}
//==============================================================================
@ -288,21 +280,13 @@ get_macro_xs(int index, int xstype, int gin, const int* gout,
{
int gout_c;
const int* gout_c_p;
int dg_c;
const int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return data::macro_xs[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
return data::macro_xs[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg);
}
//==============================================================================

6
src/state_point.cpp
View file

@ -141,13 +141,13 @@ openmc_statepoint_write(const char* filename, bool* write_source)
std::vector<int32_t> filter_ids;
filter_ids.reserve(model::tally_filters.size());
for (const auto& filt : model::tally_filters)
filter_ids.push_back(filt->id_);
filter_ids.push_back(filt->id());
write_attribute(filters_group, "ids", filter_ids);
// Write info for each filter
for (const auto& filt : model::tally_filters) {
hid_t filter_group = create_group(filters_group,
"filter " + std::to_string(filt->id_));
"filter " + std::to_string(filt->id()));
filt->to_statepoint(filter_group);
close_group(filter_group);
}
@ -188,7 +188,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
std::vector<int32_t> filter_ids;
filter_ids.reserve(tally.filters().size());
for (auto i_filt : tally.filters())
filter_ids.push_back(model::tally_filters[i_filt]->id_);
filter_ids.push_back(model::tally_filters[i_filt]->id());
write_dataset(tally_group, "filters", filter_ids);
}

1
src/summary.cpp
View file

@ -92,6 +92,7 @@ void write_geometry(hid_t file)
#ifdef DAGMC
if (settings::dagmc) {
write_attribute(geom_group, "dagmc", 1);
close_group(geom_group);
return;
}
#endif

78
src/surface.cpp
View file

@ -271,17 +271,8 @@ Direction DAGSurface::reflect(Position r, Direction u) const
return simulation::last_dir;
}
BoundingBox DAGSurface::bounding_box() const
{
moab::ErrorCode rval;
moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double min[3], max[3];
rval = dagmc_ptr_->getobb(surf, min, max);
MB_CHK_ERR_CONT(rval);
return {min[0], max[0], min[1], max[1], min[2], max[2]};
}
void DAGSurface::to_hdf5(hid_t group_id) const {}
#endif
//==============================================================================
// PeriodicSurface implementation
@ -366,9 +357,13 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
}
BoundingBox
SurfaceXPlane::bounding_box() const
SurfaceXPlane::bounding_box(bool pos_side) const
{
return {x0_, x0_, -INFTY, INFTY, -INFTY, INFTY};
if (pos_side) {
return {x0_, INFTY, -INFTY, INFTY, -INFTY, INFTY};
} else {
return {-INFTY, x0_, -INFTY, INFTY, -INFTY, INFTY};
}
}
//==============================================================================
@ -428,9 +423,13 @@ bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other,
}
BoundingBox
SurfaceYPlane::bounding_box() const
SurfaceYPlane::bounding_box(bool pos_side) const
{
return {-INFTY, INFTY, y0_, y0_, -INFTY, INFTY};
if (pos_side) {
return {-INFTY, INFTY, y0_, INFTY, -INFTY, INFTY};
} else {
return {-INFTY, INFTY, -INFTY, y0_, -INFTY, INFTY};
}
}
//==============================================================================
@ -474,9 +473,13 @@ bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other,
}
BoundingBox
SurfaceZPlane::bounding_box() const
SurfaceZPlane::bounding_box(bool pos_side) const
{
return {-INFTY, INFTY, -INFTY, INFTY, z0_, z0_};
if (pos_side) {
return {-INFTY, INFTY, -INFTY, INFTY, z0_, INFTY};
} else {
return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, z0_};
}
}
//==============================================================================
@ -538,12 +541,6 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
return false;
}
BoundingBox
SurfacePlane::bounding_box() const
{
return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
}
//==============================================================================
// Generic functions for x-, y-, and z-, cylinders
//==============================================================================
@ -645,7 +642,6 @@ Direction SurfaceXCylinder::normal(Position r) const
return axis_aligned_cylinder_normal<0, 1, 2>(r, y0_, z0_);
}
void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cylinder", false);
@ -653,6 +649,13 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const {
if (!pos_side) {
return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_};
} else {
return {};
}
}
//==============================================================================
// SurfaceYCylinder implementation
//==============================================================================
@ -686,6 +689,14 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const {
if (!pos_side) {
return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_};
} else {
return {};
}
}
//==============================================================================
// SurfaceZCylinder implementation
//==============================================================================
@ -719,6 +730,15 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const {
if (!pos_side) {
return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY};
} else {
return {};
}
}
//==============================================================================
// SurfaceSphere implementation
//==============================================================================
@ -787,6 +807,16 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
write_dataset(group_id, "coefficients", coeffs);
}
BoundingBox SurfaceSphere::bounding_box(bool pos_side) const {
if (!pos_side) {
return {x0_ - radius_, x0_ + radius_,
y0_ - radius_, y0_ + radius_,
z0_ - radius_, z0_ + radius_};
} else {
return {};
}
}
//==============================================================================
// Generic functions for x-, y-, and z-, cones
//==============================================================================
@ -1161,7 +1191,7 @@ void read_surfaces(pugi::xml_node node)
}
// See if this surface makes part of the global bounding box.
BoundingBox bb = surf->bounding_box();
auto bb = surf->bounding_box(true) & surf->bounding_box(false);
if (bb.xmin > -INFTY && bb.xmin < xmin) {
xmin = bb.xmin;
i_xmin = i_surf;

89
src/tallies/filter.cpp
View file

@ -7,6 +7,7 @@
#include "openmc/capi.h"
#include "openmc/constants.h" // for MAX_LINE_LEN;
#include "openmc/error.h"
#include "openmc/xml_interface.h"
#include "openmc/tallies/filter_azimuthal.h"
#include "openmc/tallies/filter_cell.h"
#include "openmc/tallies/filter_cellborn.h"
@ -50,8 +51,46 @@ namespace model {
// Non-member functions
//==============================================================================
Filter*
allocate_filter(const std::string& type)
extern "C" size_t tally_filters_size()
{
return model::tally_filters.size();
}
//==============================================================================
// Filter implementation
//==============================================================================
Filter::Filter() : index_{model::tally_filters.size()}
{ }
Filter::~Filter()
{
model::filter_map.erase(id_);
}
Filter* Filter::create(pugi::xml_node node)
{
// Copy filter id
if (!check_for_node(node, "id")) {
fatal_error("Must specify id for filter in tally XML file.");
}
int filter_id = std::stoi(get_node_value(node, "id"));
// Convert filter type to lower case
std::string s;
if (check_for_node(node, "type")) {
s = get_node_value(node, "type", true);
}
// Allocate according to the filter type
auto f = Filter::create(s, filter_id);
// Read filter data from XML
f->from_xml(node);
return f;
}
Filter* Filter::create(const std::string& type, int32_t id)
{
if (type == "azimuthal") {
model::tally_filters.push_back(std::make_unique<AzimuthalFilter>());
@ -100,12 +139,40 @@ allocate_filter(const std::string& type)
} else {
throw std::runtime_error{"Unknown filter type: " + type};
}
// Assign ID
model::tally_filters.back()->set_id(id);
return model::tally_filters.back().get();
}
extern "C" size_t tally_filters_size()
void Filter::set_id(int32_t id)
{
return model::tally_filters.size();
Expects(id >= -1);
// Clear entry in filter map if an ID was already assigned before
if (id_ != -1) {
model::filter_map.erase(id_);
id_ = -1;
}
// Make sure no other filter has same ID
if (model::filter_map.find(id) != model::filter_map.end()) {
throw std::runtime_error{"Two filters have the same ID: " + std::to_string(id)};
}
// If no ID specified, auto-assign next ID in sequence
if (id == -1) {
id = 0;
for (const auto& f : model::tally_filters) {
id = std::max(id, f->id_);
}
++id;
}
// Update ID and entry in filter map
id_ = id;
model::filter_map[id] = index_;
}
//==============================================================================
@ -126,7 +193,7 @@ openmc_filter_get_id(int32_t index, int32_t* id)
{
if (int err = verify_filter(index)) return err;
*id = model::tally_filters[index]->id_;
*id = model::tally_filters[index]->id();
return 0;
}
@ -135,13 +202,7 @@ openmc_filter_set_id(int32_t index, int32_t id)
{
if (int err = verify_filter(index)) return err;
if (model::filter_map.find(id) != model::filter_map.end()) {
set_errmsg("Two filters have the same ID: " + std::to_string(id));
return OPENMC_E_INVALID_ID;
}
model::tally_filters[index]->id_ = id;
model::filter_map[id] = index;
model::tally_filters[index]->set_id(id);
return 0;
}
@ -172,7 +233,7 @@ openmc_get_filter_next_id(int32_t* id)
{
int32_t largest_filter_id = 0;
for (const auto& t : model::tally_filters) {
largest_filter_id = std::max(largest_filter_id, t->id_);
largest_filter_id = std::max(largest_filter_id, t->id());
}
*id = largest_filter_id + 1;
}
@ -181,7 +242,7 @@ extern "C" int
openmc_new_filter(const char* type, int32_t* index)
{
*index = model::tally_filters.size();
allocate_filter(type);
Filter::create(type);
return 0;
}

33
src/tallies/filter_azimuthal.cpp
View file

@ -15,22 +15,35 @@ AzimuthalFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(node, "bins");
if (bins.size() > 1) {
bins_ = bins;
} else {
if (bins.size() == 1) {
// Allow a user to input a lone number which will mean that you subdivide
// [-pi,pi) evenly with the input being the number of bins
int n_angle = bins[0];
if (n_angle <= 1) fatal_error("Number of bins for azimuthal filter must "
"be greater than 1.");
if (n_angle <= 1) throw std::runtime_error{
"Number of bins for azimuthal filter must be greater than 1."};
double d_angle = 2.0 * PI / n_angle;
bins_.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins_[i] = -PI + i * d_angle;
bins_[n_angle] = PI;
bins.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins[i] = -PI + i * d_angle;
bins[n_angle] = PI;
}
this->set_bins(bins);
}
void AzimuthalFilter::set_bins(gsl::span<double> bins)
{
// Clear existing bins
bins_.clear();
bins_.reserve(bins.size());
// Copy bins, ensuring they are valid
for (gsl::index i = 0; i < bins.size(); ++i) {
if (i > 0 && bins[i] <= bins[i-1]) {
throw std::runtime_error{"Azimuthal bins must be monotonically increasing."};
}
bins_.push_back(bins[i]);
}
n_bins_ = bins_.size() - 1;

50
src/tallies/filter_cell.cpp
View file

@ -12,29 +12,39 @@ namespace openmc {
void
CellFilter::from_xml(pugi::xml_node node)
{
cells_ = get_node_array<int32_t>(node, "bins");
n_bins_ = cells_.size();
// Get cell IDs and convert into indices into the global cells vector
auto cells = get_node_array<int32_t>(node, "bins");
for (auto& c : cells) {
auto search = model::cell_map.find(c);
if (search == model::cell_map.end()) {
std::stringstream err_msg;
err_msg << "Could not find cell " << c
<< " specified on tally filter.";
throw std::runtime_error{err_msg.str()};
}
c = search->second;
}
this->set_cells(cells);
}
void
CellFilter::initialize()
CellFilter::set_cells(gsl::span<int32_t> cells)
{
// Convert cell IDs to indices of the global array.
for (auto& c : cells_) {
auto search = model::cell_map.find(c);
if (search != model::cell_map.end()) {
c = search->second;
} else {
std::stringstream err_msg;
err_msg << "Could not find cell " << c << " specified on tally filter.";
fatal_error(err_msg);
}
// Clear existing cells
cells_.clear();
cells_.reserve(cells.size());
map_.clear();
// Update cells and mapping
for (auto& index : cells) {
Expects(index >= 0);
Expects(index < model::cells.size());
cells_.push_back(index);
map_[index] = cells_.size() - 1;
}
// Populate the index->bin map.
for (int i = 0; i < cells_.size(); i++) {
map_[cells_[i]] = i;
}
n_bins_ = cells_.size();
}
void
@ -70,7 +80,7 @@ CellFilter::text_label(int bin) const
//==============================================================================
extern "C" int
openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n)
openmc_cell_filter_get_bins(int32_t index, const int32_t** cells, int32_t* n)
{
if (int err = verify_filter(index)) return err;
@ -81,8 +91,8 @@ openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n)
}
auto cell_filt = static_cast<CellFilter*>(filt);
*cells = cell_filt->cells_.data();
*n = cell_filt->cells_.size();
*cells = cell_filt->cells().data();
*n = cell_filt->cells().size();
return 0;
}

25
src/tallies/filter_delayedgroup.cpp
View file

@ -8,21 +8,32 @@ namespace openmc {
void
DelayedGroupFilter::from_xml(pugi::xml_node node)
{
groups_ = get_node_array<int>(node, "bins");
n_bins_ = groups_.size();
auto groups = get_node_array<int>(node, "bins");
this->set_groups(groups);
}
void
DelayedGroupFilter::set_groups(gsl::span<int> groups)
{
// Clear existing groups
groups_.clear();
groups_.reserve(groups.size());
// Make sure all the group index values are valid.
// TODO: do these need to be decremented for zero-based indexing?
for (auto group : groups_) {
for (auto group : groups) {
if (group < 1) {
fatal_error("Encountered delayedgroup bin with index "
+ std::to_string(group) + " which is less than 1");
throw std::invalid_argument{"Encountered delayedgroup bin with index "
+ std::to_string(group) + " which is less than 1"};
} else if (group > MAX_DELAYED_GROUPS) {
fatal_error("Encountered delayedgroup bin with index "
throw std::invalid_argument{"Encountered delayedgroup bin with index "
+ std::to_string(group) + " which is greater than MAX_DELATED_GROUPS ("
+ std::to_string(MAX_DELAYED_GROUPS) + ")");
+ std::to_string(MAX_DELAYED_GROUPS) + ")"};
}
groups_.push_back(group);
}
n_bins_ = groups_.size();
}
void

27
src/tallies/filter_distribcell.cpp
View file

@ -15,23 +15,26 @@ DistribcellFilter::from_xml(pugi::xml_node node)
if (cells.size() != 1) {
fatal_error("Only one cell can be specified per distribcell filter.");
}
cell_ = cells[0];
}
void
DistribcellFilter::initialize()
{
// Convert the cell ID to an index of the global array.
auto search = model::cell_map.find(cell_);
if (search != model::cell_map.end()) {
cell_ = search->second;
n_bins_ = model::cells[cell_]->n_instances_;
} else {
// Find index in global cells vector corresponding to cell ID
auto search = model::cell_map.find(cells[0]);
if (search == model::cell_map.end()) {
std::stringstream err_msg;
err_msg << "Could not find cell " << cell_
<< " specified on tally filter.";
fatal_error(err_msg);
throw std::runtime_error{err_msg.str()};
}
this->set_cell(search->second);
}
void
DistribcellFilter::set_cell(int32_t cell)
{
Expects(cell >= 0);
Expects(cell < model::cells.size());
cell_ = cell;
n_bins_ = model::cells[cell]->n_instances_;
}
void

35
src/tallies/filter_energy.cpp
View file

@ -16,7 +16,25 @@ namespace openmc {
void
EnergyFilter::from_xml(pugi::xml_node node)
{
bins_ = get_node_array<double>(node, "bins");
auto bins = get_node_array<double>(node, "bins");
this->set_bins(bins);
}
void
EnergyFilter::set_bins(gsl::span<const double> bins)
{
// Clear existing bins
bins_.clear();
bins_.reserve(bins.size());
// Copy bins, ensuring they are valid
for (gsl::index i = 0; i < bins.size(); ++i) {
if (i > 0 && bins[i] <= bins[i-1]) {
throw std::runtime_error{"Energy bins must be monotonically increasing."};
}
bins_.push_back(bins[i]);
}
n_bins_ = bins_.size() - 1;
// In MG mode, check if the filter bins match the transport bins.
@ -27,7 +45,7 @@ EnergyFilter::from_xml(pugi::xml_node node)
if (!settings::run_CE) {
if (n_bins_ == data::num_energy_groups) {
matches_transport_groups_ = true;
for (auto i = 0; i < n_bins_ + 1; i++) {
for (gsl::index i = 0; i < n_bins_ + 1; ++i) {
if (data::rev_energy_bins[i] != bins_[i]) {
matches_transport_groups_ = false;
break;
@ -111,7 +129,7 @@ EnergyoutFilter::text_label(int bin) const
//==============================================================================
extern"C" int
openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n)
openmc_energy_filter_get_bins(int32_t index, const double** energies, size_t* n)
{
// Make sure this is a valid index to an allocated filter.
if (int err = verify_filter(index)) return err;
@ -127,13 +145,13 @@ openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n)
}
// Output the bins.
*energies = filt->bins_.data();
*n = filt->bins_.size();
*energies = filt->bins().data();
*n = filt->bins().size();
return 0;
}
extern "C" int
openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
openmc_energy_filter_set_bins(int32_t index, size_t n, const double* energies)
{
// Make sure this is a valid index to an allocated filter.
if (int err = verify_filter(index)) return err;
@ -149,10 +167,7 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
}
// Update the filter.
filt->bins_.clear();
filt->bins_.resize(n);
for (int i = 0; i < n; i++) filt->bins_[i] = energies[i];
filt->n_bins_ = n - 1;
filt->set_bins({energies, n});
return 0;
}

97
src/tallies/filter_energyfunc.cpp
View file

@ -21,12 +21,37 @@ EnergyFunctionFilter::from_xml(pugi::xml_node node)
if (!check_for_node(node, "energy"))
fatal_error("Energy grid not specified for EnergyFunction filter.");
energy_ = get_node_array<double>(node, "energy");
auto energy = get_node_array<double>(node, "energy");
if (!check_for_node(node, "y"))
fatal_error("y values not specified for EnergyFunction filter.");
y_ = get_node_array<double>(node, "y");
auto y = get_node_array<double>(node, "y");
this->set_data(energy, y);
}
void
EnergyFunctionFilter::set_data(gsl::span<const double> energy,
gsl::span<const double> y)
{
// Check for consistent sizes with new data
if (energy.size() != y.size()) {
fatal_error("Energy grid and y values are not consistent");
}
energy_.clear();
energy_.reserve(energy.size());
y_.clear();
y_.reserve(y.size());
// Copy over energy values, ensuring they are valid
for (gsl::index i = 0; i < energy.size(); ++i) {
if (i > 0 && energy[i] <= energy[i - 1]) {
throw std::runtime_error{"Energy bins must be monotonically increasing."};
}
energy_.push_back(energy[i]);
y_.push_back(y[i]);
}
}
void
@ -65,4 +90,72 @@ EnergyFunctionFilter::text_label(int bin) const
return out.str();
}
//==============================================================================
// C-API functions
//==============================================================================
extern "C" int
openmc_energyfunc_filter_set_data(int32_t index, size_t n, const double* energy,
const double* y)
{
// Ensure this is a valid index to allocated filter
if (int err = verify_filter(index)) return err;
// Get a pointer to the filter
const auto& filt_base = model::tally_filters[index].get();
// Downcast to EnergyFunctionFilter
auto* filt = dynamic_cast<EnergyFunctionFilter*>(filt_base);
// Check if a valid filter was produced
if (!filt) {
set_errmsg("Tried to set interpolation data for non-energy function filter.");
return OPENMC_E_INVALID_TYPE;
}
filt->set_data({energy, n}, {y, n});
return 0;
}
extern "C" int
openmc_energyfunc_filter_get_energy(int32_t index, size_t *n, const double** energy)
{
// ensure this is a valid index to allocated filter
if (int err = verify_filter(index)) return err;
// get a pointer to the filter
const auto& filt_base = model::tally_filters[index].get();
// downcast to EnergyFunctionFilter
auto* filt = dynamic_cast<EnergyFunctionFilter*>(filt_base);
// check if a valid filter was produced
if (!filt) {
set_errmsg("Tried to set interpolation data for non-energy function filter.");
return OPENMC_E_INVALID_TYPE;
}
*energy = filt->energy().data();
*n = filt->energy().size();
return 0;
}
extern "C" int
openmc_energyfunc_filter_get_y(int32_t index, size_t *n, const double** y)
{
// ensure this is a valid index to allocated filter
if (int err = verify_filter(index)) return err;
// get a pointer to the filter
const auto& filt_base = model::tally_filters[index].get();
// downcast to EnergyFunctionFilter
auto* filt = dynamic_cast<EnergyFunctionFilter*>(filt_base);
// check if a valid filter was produced
if (!filt) {
set_errmsg("Tried to set interpolation data for non-energy function filter.");
return OPENMC_E_INVALID_TYPE;
}
*y = filt->y().data();
*n = filt->y().size();
return 0;
}
} // namespace openmc

16
src/tallies/filter_legendre.cpp
View file

@ -10,7 +10,16 @@ namespace openmc {
void
LegendreFilter::from_xml(pugi::xml_node node)
{
order_ = std::stoi(get_node_value(node, "order"));
this->set_order(std::stoi(get_node_value(node, "order")));
}
void
LegendreFilter::set_order(int order)
{
if (order < 0) {
throw std::invalid_argument{"Legendre order must be non-negative."};
}
order_ = order;
n_bins_ = order_ + 1;
}
@ -60,7 +69,7 @@ openmc_legendre_filter_get_order(int32_t index, int* order)
}
// Output the order.
*order = filt->order_;
*order = filt->order();
return 0;
}
@ -81,8 +90,7 @@ openmc_legendre_filter_set_order(int32_t index, int order)
}
// Update the filter.
filt->order_ = order;
filt->n_bins_ = order + 1;
filt->set_order(order);
return 0;
}

57
src/tallies/filter_material.cpp
View file

@ -3,7 +3,6 @@
#include <sstream>
#include "openmc/capi.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/xml_interface.h"
@ -12,30 +11,39 @@ namespace openmc {
void
MaterialFilter::from_xml(pugi::xml_node node)
{
materials_ = get_node_array<int32_t>(node, "bins");
n_bins_ = materials_.size();
}
void
MaterialFilter::initialize()
{
// Convert material IDs to indices of the global array.
for (auto& m : materials_) {
// Get material IDs and convert to indices in the global materials vector
auto mats = get_node_array<int32_t>(node, "bins");
for (auto& m : mats) {
auto search = model::material_map.find(m);
if (search != model::material_map.end()) {
m = search->second;
} else {
if (search == model::material_map.end()) {
std::stringstream err_msg;
err_msg << "Could not find material " << m
<< " specified on tally filter.";
fatal_error(err_msg);
throw std::runtime_error{err_msg.str()};
}
m = search->second;
}
// Populate the index->bin map.
for (int i = 0; i < materials_.size(); i++) {
map_[materials_[i]] = i;
this->set_materials(mats);
}
void
MaterialFilter::set_materials(gsl::span<const int32_t> materials)
{
// Clear existing materials
materials_.clear();
materials_.reserve(materials.size());
map_.clear();
// Update materials and mapping
for (auto& index : materials) {
Expects(index >= 0);
Expects(index < model::materials.size());
materials_.push_back(index);
map_[index] = materials_.size() - 1;
}
n_bins_ = materials_.size();
}
void
@ -69,7 +77,7 @@ MaterialFilter::text_label(int bin) const
//==============================================================================
extern "C" int
openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n)
openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n)
{
// Make sure this is a valid index to an allocated filter.
if (int err = verify_filter(index)) return err;
@ -85,13 +93,13 @@ openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n)
}
// Output the bins.
*bins = filt->materials_.data();
*n = filt->materials_.size();
*bins = filt->materials().data();
*n = filt->materials().size();
return 0;
}
extern "C" int
openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins)
openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins)
{
// Make sure this is a valid index to an allocated filter.
if (int err = verify_filter(index)) return err;
@ -107,12 +115,7 @@ openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins)
}
// Update the filter.
filt->materials_.clear();
filt->materials_.resize(n);
for (int i = 0; i < n; i++) filt->materials_[i] = bins[i];
filt->n_bins_ = filt->materials_.size();
filt->map_.clear();
for (int i = 0; i < n; i++) filt->map_[filt->materials_[i]] = i;
filt->set_materials({bins, n});
return 0;
}

34
src/tallies/filter_mu.cpp
View file

@ -13,22 +13,36 @@ MuFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(node, "bins");
if (bins.size() > 1) {
bins_ = bins;
} else {
if (bins.size() == 1) {
// Allow a user to input a lone number which will mean that you subdivide
// [-1,1) evenly with the input being the number of bins
int n_angle = bins[0];
if (n_angle <= 1) fatal_error("Number of bins for mu filter must "
"be greater than 1.");
if (n_angle <= 1) throw std::runtime_error{
"Number of bins for mu filter must be greater than 1."};
double d_angle = 2.0 / n_angle;
bins_.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins_[i] = -1 + i * d_angle;
bins_[n_angle] = 1;
bins.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins[i] = -1 + i * d_angle;
bins[n_angle] = 1;
}
this->set_bins(bins);
}
void
MuFilter::set_bins(gsl::span<double> bins)
{
// Clear existing bins
bins_.clear();
bins_.reserve(bins.size());
// Copy bins, ensuring they are valid
for (gsl::index i = 0; i < bins.size(); ++i) {
if (i > 0 && bins[i] <= bins[i-1]) {
throw std::runtime_error{"Mu bins must be monotonically increasing."};
}
bins_.push_back(bins[i]);
}
n_bins_ = bins_.size() - 1;

46
src/tallies/filter_particle.cpp
View file

@ -7,9 +7,34 @@ namespace openmc {
void
ParticleFilter::from_xml(pugi::xml_node node)
{
auto particles = get_node_array<int>(node, "bins");
auto particles = get_node_array<std::string>(node, "bins");
// Convert to vector of Particle::Type
std::vector<Particle::Type> types;
for (auto& p : particles) {
particles_.push_back(static_cast<Particle::Type>(p - 1));
if (p == "neutron") {
types.push_back(Particle::Type::neutron);
} else if (p == "photon") {
types.push_back(Particle::Type::photon);
} else if (p == "electron") {
types.push_back(Particle::Type::electron);
} else if (p == "positron") {
types.push_back(Particle::Type::positron);
}
}
this->set_particles(types);
}
void
ParticleFilter::set_particles(gsl::span<Particle::Type> particles)
{
// Clear existing particles
particles_.clear();
particles_.reserve(particles.size());
// Set particles and number of bins
for (auto p : particles) {
particles_.push_back(p);
}
n_bins_ = particles_.size();
}
@ -30,9 +55,22 @@ void
ParticleFilter::to_statepoint(hid_t filter_group) const
{
Filter::to_statepoint(filter_group);
std::vector<int> particles;
std::vector<std::string> particles;
for (auto p : particles_) {
particles.push_back(static_cast<int>(p) + 1);
switch (p) {
case Particle::Type::neutron:
particles.push_back("neutron");
break;
case Particle::Type::photon:
particles.push_back("photon");
break;
case Particle::Type::electron:
particles.push_back("electron");
break;
case Particle::Type::positron:
particles.push_back("positron");
break;
}
}
write_dataset(filter_group, "bins", particles);
}

34
src/tallies/filter_polar.cpp
View file

@ -14,22 +14,36 @@ PolarFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(node, "bins");
if (bins.size() > 1) {
bins_ = bins;
} else {
if (bins.size() == 1) {
// Allow a user to input a lone number which will mean that you subdivide
// [0,pi] evenly with the input being the number of bins
int n_angle = bins[0];
if (n_angle <= 1) fatal_error("Number of bins for polar filter must "
"be greater than 1.");
if (n_angle <= 1) throw std::runtime_error{
"Number of bins for polar filter must be greater than 1."};
double d_angle = PI / n_angle;
bins_.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins_[i] = i * d_angle;
bins_[n_angle] = PI;
bins.resize(n_angle + 1);
for (int i = 0; i < n_angle; i++) bins[i] = i * d_angle;
bins[n_angle] = PI;
}
this->set_bins(bins);
}
void
PolarFilter::set_bins(gsl::span<double> bins)
{
// Clear existing bins
bins_.clear();
bins_.reserve(bins.size());
// Copy bins, ensuring they are valid
for (gsl::index i = 0; i < bins.size(); ++i) {
if (i > 0 && bins[i] <= bins[i-1]) {
throw std::runtime_error{"Polar bins must be monotonically increasing."};
}
bins_.push_back(bins[i]);
}
n_bins_ = bins_.size() - 1;

58
src/tallies/filter_sph_harm.cpp
View file

@ -12,21 +12,34 @@ namespace openmc {
void
SphericalHarmonicsFilter::from_xml(pugi::xml_node node)
{
order_ = std::stoi(get_node_value(node, "order"));
n_bins_ = (order_ + 1) * (order_ + 1);
this->set_order(std::stoi(get_node_value(node, "order")));
if (check_for_node(node, "cosine")) {
auto cos = get_node_value(node, "cosine", true);
if (cos == "scatter") {
cosine_ = SphericalHarmonicsCosine::scatter;
} else if (cos == "particle") {
cosine_ = SphericalHarmonicsCosine::particle;
} else {
std::stringstream err_msg;
err_msg << "Unrecognized cosine type, \"" << cos
<< "\" in spherical harmonics filter";
fatal_error(err_msg);
}
this->set_cosine(get_node_value(node, "cosine", true));
}
}
void
SphericalHarmonicsFilter::set_order(int order)
{
if (order < 0) {
throw std::invalid_argument{"Spherical harmonics order must be non-negative."};
}
order_ = order;
n_bins_ = (order_ + 1) * (order_ + 1);
}
void
SphericalHarmonicsFilter::set_cosine(gsl::cstring_span cosine)
{
if (cosine == "scatter") {
cosine_ = SphericalHarmonicsCosine::scatter;
} else if (cosine == "particle") {
cosine_ = SphericalHarmonicsCosine::particle;
} else {
std::stringstream err_msg;
err_msg << "Unrecognized cosine type, \"" << cosine
<< "\" in spherical harmonics filter";
throw std::invalid_argument{err_msg.str()};
}
}
@ -121,7 +134,7 @@ openmc_sphharm_filter_get_order(int32_t index, int* order)
if (err) return err;
// Output the order.
*order = filt->order_;
*order = filt->order();
return 0;
}
@ -135,7 +148,7 @@ openmc_sphharm_filter_get_cosine(int32_t index, char cosine[])
if (err) return err;
// Output the cosine.
if (filt->cosine_ == SphericalHarmonicsCosine::scatter) {
if (filt->cosine() == SphericalHarmonicsCosine::scatter) {
strcpy(cosine, "scatter");
} else {
strcpy(cosine, "particle");
@ -153,8 +166,7 @@ openmc_sphharm_filter_set_order(int32_t index, int order)
if (err) return err;
// Update the filter.
filt->order_ = order;
filt->n_bins_ = (order + 1) * (order + 1);
filt->set_order(order);
return 0;
}
@ -168,12 +180,10 @@ openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[])
if (err) return err;
// Update the filter.
if (strcmp(cosine, "scatter") == 0) {
filt->cosine_ = SphericalHarmonicsCosine::scatter;
} else if (strcmp(cosine, "particle") == 0) {
filt->cosine_ = SphericalHarmonicsCosine::particle;
} else {
set_errmsg("Invalid spherical harmonics cosine.");
try {
filt->set_cosine(cosine);
} catch (const std::invalid_argument& e) {
set_errmsg(e.what());
return OPENMC_E_INVALID_ARGUMENT;
}
return 0;

67
src/tallies/filter_sptl_legendre.cpp
View file

@ -12,25 +12,54 @@ namespace openmc {
void
SpatialLegendreFilter::from_xml(pugi::xml_node node)
{
order_ = std::stoi(get_node_value(node, "order"));
this->set_order(std::stoi(get_node_value(node, "order")));
auto axis = get_node_value(node, "axis");
if (axis == "x") {
axis_ = LegendreAxis::x;
} else if (axis == "y") {
axis_ = LegendreAxis::y;
} else if (axis == "z") {
axis_ = LegendreAxis::z;
} else {
fatal_error("Unrecognized axis on SpatialLegendreFilter");
switch (axis[0]) {
case 'x':
this->set_axis(LegendreAxis::x);
break;
case 'y':
this->set_axis(LegendreAxis::y);
break;
case 'z':
this->set_axis(LegendreAxis::z);
break;
default:
throw std::runtime_error{"Axis for SpatialLegendreFilter must be 'x', 'y', or 'z'"};
}
min_ = std::stod(get_node_value(node, "min"));
max_ = std::stod(get_node_value(node, "max"));
double min = std::stod(get_node_value(node, "min"));
double max = std::stod(get_node_value(node, "max"));
this->set_minmax(min, max);
}
void
SpatialLegendreFilter::set_order(int order)
{
if (order < 0) {
throw std::invalid_argument{"Legendre order must be non-negative."};
}
order_ = order;
n_bins_ = order_ + 1;
}
void
SpatialLegendreFilter::set_axis(LegendreAxis axis)
{
axis_ = axis;
}
void
SpatialLegendreFilter::set_minmax(double min, double max)
{
if (max < min) {
throw std::invalid_argument{"Maximum value must be greater than minimum value"};
}
min_ = min;
max_ = max;
}
void
SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
@ -126,7 +155,7 @@ openmc_spatial_legendre_filter_get_order(int32_t index, int* order)
if (err) return err;
// Output the order.
*order = filt->order_;
*order = filt->order();
return 0;
}
@ -141,9 +170,9 @@ openmc_spatial_legendre_filter_get_params(int32_t index, int* axis,
if (err) return err;
// Output the params.
*axis = static_cast<int>(filt->axis_);
*min = filt->min_;
*max = filt->max_;
*axis = static_cast<int>(filt->axis());
*min = filt->min();
*max = filt->max();
return 0;
}
@ -157,8 +186,7 @@ openmc_spatial_legendre_filter_set_order(int32_t index, int order)
if (err) return err;
// Update the filter.
filt->order_ = order;
filt->n_bins_ = order + 1;
filt->set_order(order);
return 0;
}
@ -173,9 +201,8 @@ openmc_spatial_legendre_filter_set_params(int32_t index, const int* axis,
if (err) return err;
// Update the filter.
if (axis) filt->axis_ = static_cast<LegendreAxis>(*axis);
if (min) filt->min_ = *min;
if (max) filt->max_ = *max;
if (axis) filt->set_axis(static_cast<LegendreAxis>(*axis));
if (min && max) filt->set_minmax(*min, *max);
return 0;
}

41
src/tallies/filter_surface.cpp
View file

@ -11,30 +11,41 @@ namespace openmc {
void
SurfaceFilter::from_xml(pugi::xml_node node)
{
surfaces_ = get_node_array<int32_t>(node, "bins");
n_bins_ = surfaces_.size();
}
auto surfaces = get_node_array<int32_t>(node, "bins");
void
SurfaceFilter::initialize()
{
// Convert surface IDs to indices of the global array.
for (auto& s : surfaces_) {
// Convert surface IDs to indices of the global surfaces vector.
for (auto& s : surfaces) {
auto search = model::surface_map.find(s);
if (search != model::surface_map.end()) {
s = search->second;
} else {
if (search == model::surface_map.end()) {
std::stringstream err_msg;
err_msg << "Could not find surface " << s
<< " specified on tally filter.";
fatal_error(err_msg);
throw std::runtime_error{err_msg.str()};
}
s = search->second;
}
// Populate the index->bin map.
for (int i = 0; i < surfaces_.size(); i++) {
map_[surfaces_[i]] = i;
this->set_surfaces(surfaces);
}
void
SurfaceFilter::set_surfaces(gsl::span<int32_t> surfaces)
{
// Clear existing surfaces
surfaces_.clear();
surfaces_.reserve(surfaces.size());
map_.clear();
// Update surfaces and mapping
for (auto& index : surfaces) {
Expects(index >= 0);
Expects(index < model::surfaces.size());
surfaces_.push_back(index);
map_[index] = surfaces_.size() - 1;
}
n_bins_ = surfaces_.size();
}
void

41
src/tallies/filter_universe.cpp
View file

@ -11,30 +11,39 @@ namespace openmc {
void
UniverseFilter::from_xml(pugi::xml_node node)
{
universes_ = get_node_array<int32_t>(node, "bins");
n_bins_ = universes_.size();
}
void
UniverseFilter::initialize()
{
// Convert universe IDs to indices of the global array.
for (auto& u : universes_) {
// Get material IDs and convert to indices in the global materials vector
auto universes = get_node_array<int32_t>(node, "bins");
for (auto& u : universes) {
auto search = model::universe_map.find(u);
if (search != model::universe_map.end()) {
u = search->second;
} else {
if (search == model::universe_map.end()) {
std::stringstream err_msg;
err_msg << "Could not find universe " << u
<< " specified on tally filter.";
fatal_error(err_msg);
throw std::runtime_error{err_msg.str()};
}
u = search->second;
}
// Populate the index->bin map.
for (int i = 0; i < universes_.size(); i++) {
map_[universes_[i]] = i;
this->set_universes(universes);
}
void
UniverseFilter::set_universes(gsl::span<int32_t> universes)
{
// Clear existing universes
universes_.clear();
universes_.reserve(universes.size());
map_.clear();
// Update universes and mapping
for (auto& index : universes) {
Expects(index >= 0);
Expects(index < model::universes.size());
universes_.push_back(index);
map_[index] = universes_.size() - 1;
}
n_bins_ = universes_.size();
}
void

17
src/tallies/filter_zernike.cpp
View file

@ -74,6 +74,9 @@ ZernikeFilter::text_label(int bin) const
void
ZernikeFilter::set_order(int order)
{
if (order < 0) {
throw std::invalid_argument{"Zernike order must be non-negative."};
}
order_ = order;
n_bins_ = ((order+1) * (order+2)) / 2;
}
@ -111,7 +114,7 @@ ZernikeRadialFilter::text_label(int bin) const
void
ZernikeRadialFilter::set_order(int order)
{
order_ = order;
ZernikeFilter::set_order(order);
n_bins_ = order / 2 + 1;
}
@ -165,9 +168,9 @@ openmc_zernike_filter_get_params(int32_t index, double* x, double* y,
if (err) return err;
// Output the params.
*x = filt->x_;
*y = filt->y_;
*r = filt->r_;
*x = filt->x();
*y = filt->y();
*r = filt->r();
return 0;
}
@ -196,9 +199,9 @@ openmc_zernike_filter_set_params(int32_t index, const double* x,
if (err) return err;
// Update the filter.
if (x) filt->x_ = *x;
if (y) filt->y_ = *y;
if (r) filt->r_ = *r;
if (x) filt->set_x(*x);
if (y) filt->set_y(*y);
if (r) filt->set_r(*r);
return 0;
}

594
src/tallies/tally.cpp
View file

@ -240,37 +240,293 @@ score_str_to_int(std::string score_str)
// Tally object implementation
//==============================================================================
Tally::Tally()
Tally::Tally(int32_t id)
: index_{model::tallies.size()}
{
this->set_filters(nullptr, 0);
this->set_id(id);
this->set_filters({});
}
void
Tally::init_from_xml(pugi::xml_node node)
Tally::Tally(pugi::xml_node node)
: index_{model::tallies.size()}
{
// Copy and set tally id
if (!check_for_node(node, "id")) {
throw std::runtime_error{"Must specify id for tally in tally XML file."};
}
int32_t id = std::stoi(get_node_value(node, "id"));
this->set_id(id);
if (check_for_node(node, "name")) name_ = get_node_value(node, "name");
// =======================================================================
// READ DATA FOR FILTERS
// Check if user is using old XML format and throw an error if so
if (check_for_node(node, "filter")) {
throw std::runtime_error{"Tally filters must be specified independently of "
"tallies in a <filter> element. The <tally> element itself should "
"have a list of filters that apply, e.g., <filters>1 2</filters> "
"where 1 and 2 are the IDs of filters specified outside of "
"<tally>."};
}
// Determine number of filters
std::vector<int> filter_ids;
if (check_for_node(node, "filters")) {
filter_ids = get_node_array<int>(node, "filters");
}
// Allocate and store filter user ids
std::vector<Filter*> filters;
for (int filter_id : filter_ids) {
// Determine if filter ID is valid
auto it = model::filter_map.find(filter_id);
if (it == model::filter_map.end()) {
throw std::runtime_error{"Could not find filter " + std::to_string(filter_id)
+ " specified on tally " + std::to_string(id_)};
}
// Store the index of the filter
filters.push_back(model::tally_filters[it->second].get());
}
// Set the filters
this->set_filters(filters);
// Check for the presence of certain filter types
bool has_energyout = energyout_filter_ >= 0;
int particle_filter_index = C_NONE;
for (gsl::index j = 0; j < filters_.size(); ++j) {
int i_filter = filters_[j];
const auto& f = model::tally_filters[i_filter].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
if (pf) particle_filter_index = i_filter;
// Change the tally estimator if a filter demands it
std::string filt_type = f->type();
if (filt_type == "energyout" || filt_type == "legendre") {
estimator_ = ESTIMATOR_ANALOG;
} else if (filt_type == "sphericalharmonics") {
auto sf = dynamic_cast<SphericalHarmonicsFilter*>(f);
if (sf->cosine() == SphericalHarmonicsCosine::scatter) {
estimator_ = ESTIMATOR_ANALOG;
}
} else if (filt_type == "spatiallegendre" || filt_type == "zernike"
|| filt_type == "zernikeradial") {
estimator_ = ESTIMATOR_COLLISION;
}
}
// =======================================================================
// READ DATA FOR NUCLIDES
this->set_nuclides(node);
// =======================================================================
// READ DATA FOR SCORES
this->set_scores(node);
if (!check_for_node(node, "scores")) {
fatal_error("No scores specified on tally " + std::to_string(id_)
+ ".");
}
// Check if tally is compatible with particle type
if (settings::photon_transport) {
if (particle_filter_index == C_NONE) {
for (int score : scores_) {
switch (score) {
case SCORE_INVERSE_VELOCITY:
fatal_error("Particle filter must be used with photon "
"transport on and inverse velocity score");
break;
case SCORE_FLUX:
case SCORE_TOTAL:
case SCORE_SCATTER:
case SCORE_NU_SCATTER:
case SCORE_ABSORPTION:
case SCORE_FISSION:
case SCORE_NU_FISSION:
case SCORE_CURRENT:
case SCORE_EVENTS:
case SCORE_DELAYED_NU_FISSION:
case SCORE_PROMPT_NU_FISSION:
case SCORE_DECAY_RATE:
warning("Particle filter is not used with photon transport"
" on and " + reaction_name(score) + " score.");
break;
}
}
} else {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (auto p : pf->particles()) {
if (p == Particle::Type::electron ||
p == Particle::Type::positron) {
estimator_ = ESTIMATOR_ANALOG;
}
}
}
} else {
if (particle_filter_index >= 0) {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (auto p : pf->particles()) {
if (p != Particle::Type::neutron) {
warning("Particle filter other than NEUTRON used with photon "
"transport turned off. All tallies for particle type " +
std::to_string(static_cast<int>(p)) + " will have no scores");
}
}
}
}
// Check for a tally derivative.
if (check_for_node(node, "derivative")) {
int deriv_id = std::stoi(get_node_value(node, "derivative"));
// Find the derivative with the given id, and store it's index.
auto it = model::tally_deriv_map.find(deriv_id);
if (it == model::tally_deriv_map.end()) {
fatal_error("Could not find derivative " + std::to_string(deriv_id)
+ " specified on tally " + std::to_string(id_));
}
deriv_ = it->second;
// Only analog or collision estimators are supported for differential
// tallies.
if (estimator_ == ESTIMATOR_TRACKLENGTH) {
estimator_ = ESTIMATOR_COLLISION;
}
const auto& deriv = model::tally_derivs[deriv_];
if (deriv.variable == DIFF_NUCLIDE_DENSITY
|| deriv.variable == DIFF_TEMPERATURE) {
for (int i_nuc : nuclides_) {
if (has_energyout && i_nuc == -1) {
fatal_error("Error on tally " + std::to_string(id_)
+ ": Cannot use a 'nuclide_density' or 'temperature' "
"derivative on a tally with an outgoing energy filter and "
"'total' nuclide rate. Instead, tally each nuclide in the "
"material individually.");
// Note that diff tallies with these characteristics would work
// correctly if no tally events occur in the perturbed material
// (e.g. pertrubing moderator but only tallying fuel), but this
// case would be hard to check for by only reading inputs.
}
}
}
}
// If settings.xml trigger is turned on, create tally triggers
if (settings::trigger_on) {
this->init_triggers(node);
}
// =======================================================================
// SET TALLY ESTIMATOR
// Check if user specified estimator
if (check_for_node(node, "estimator")) {
std::string est = get_node_value(node, "estimator");
if (est == "analog") {
estimator_ = ESTIMATOR_ANALOG;
} else if (est == "tracklength" || est == "track-length"
|| est == "pathlength" || est == "path-length") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (estimator_ == ESTIMATOR_ANALOG) {
throw std::runtime_error{"Cannot use track-length estimator for tally "
+ std::to_string(id_)};
}
// Set estimator to track-length estimator
estimator_ = ESTIMATOR_TRACKLENGTH;
} else if (est == "collision") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (estimator_ == ESTIMATOR_ANALOG) {
throw std::runtime_error{"Cannot use collision estimator for tally " +
std::to_string(id_)};
}
// Set estimator to collision estimator
estimator_ = ESTIMATOR_COLLISION;
} else {
throw std::runtime_error{"Invalid estimator '" + est + "' on tally " +
std::to_string(id_)};
}
}
}
Tally::~Tally()
{
model::tally_map.erase(id_);
}
Tally*
Tally::create(int32_t id)
{
model::tallies.push_back(std::make_unique<Tally>(id));
return model::tallies.back().get();
}
void
Tally::set_filters(const int32_t filter_indices[], int n)
Tally::set_id(int32_t id)
{
Expects(id >= -1);
// Clear entry in tally map if an ID was already assigned before
if (id_ != -1) {
model::tally_map.erase(id_);
id_ = -1;
}
// Make sure no other tally has the same ID
if (model::tally_map.find(id) != model::tally_map.end()) {
throw std::runtime_error{"Two tallies have the same ID: " + std::to_string(id)};
}
// If no ID specified, auto-assign next ID in sequence
if (id == -1) {
id = 0;
for (const auto& t : model::tallies) {
id = std::max(id, t->id_);
}
++id;
}
// Update ID and entry in tally map
id_ = id;
model::tally_map[id] = index_;
}
void
Tally::set_filters(gsl::span<Filter*> filters)
{
// Clear old data.
filters_.clear();
strides_.clear();
// Copy in the given filter indices.
filters_.assign(filter_indices, filter_indices + n);
auto n = filters.size();
filters_.reserve(n);
for (int i = 0; i < n; ++i) {
auto i_filt = filters_[i];
if (i_filt < 0 || i_filt >= model::tally_filters.size())
throw std::out_of_range("Index in tally filter array out of bounds.");
// Add index to vector of filters
auto& f {filters[i]};
filters_.push_back(model::filter_map.at(f->id()));
// Keep track of indices for special filters.
const auto* filt = model::tally_filters[i_filt].get();
if (dynamic_cast<const EnergyoutFilter*>(filt)) {
if (dynamic_cast<const EnergyoutFilter*>(f)) {
energyout_filter_ = i;
} else if (dynamic_cast<const DelayedGroupFilter*>(filt)) {
} else if (dynamic_cast<const DelayedGroupFilter*>(f)) {
delayedgroup_filter_ = i;
}
}
@ -282,7 +538,7 @@ Tally::set_filters(const int32_t filter_indices[], int n)
int stride = 1;
for (int i = n-1; i >= 0; --i) {
strides_[i] = stride;
stride *= model::tally_filters[filters_[i]]->n_bins_;
stride *= model::tally_filters[filters_[i]]->n_bins();
}
n_filter_bins_ = stride;
}
@ -298,7 +554,7 @@ Tally::set_scores(pugi::xml_node node)
}
void
Tally::set_scores(std::vector<std::string> scores)
Tally::set_scores(const std::vector<std::string>& scores)
{
// Reset state and prepare for the new scores.
scores_.clear();
@ -450,17 +706,25 @@ Tally::set_nuclides(pugi::xml_node node)
// The user provided specifics nuclides. Parse it as an array with either
// "total" or a nuclide name like "U-235" in each position.
auto words = get_node_array<std::string>(node, "nuclides");
for (auto word : words) {
if (word == "total") {
nuclides_.push_back(-1);
} else {
auto search = data::nuclide_map.find(word);
if (search == data::nuclide_map.end())
fatal_error("Could not find the nuclide " + word
+ " specified in tally " + std::to_string(id_)
+ " in any material");
nuclides_.push_back(search->second);
}
this->set_nuclides(words);
}
}
void
Tally::set_nuclides(const std::vector<std::string>& nuclides)
{
nuclides_.clear();
for (const auto& nuc : nuclides) {
if (nuc == "total") {
nuclides_.push_back(-1);
} else {
auto search = data::nuclide_map.find(nuc);
if (search == data::nuclide_map.end())
fatal_error("Could not find the nuclide " + nuc
+ " specified in tally " + std::to_string(id_)
+ " in any material");
nuclides_.push_back(search->second);
}
}
}
@ -614,36 +878,7 @@ void read_tallies_xml()
// Check for user filters and allocate
for (auto node_filt : root.children("filter")) {
// Copy filter id
if (!check_for_node(node_filt, "id")) {
fatal_error("Must specify id for filter in tally XML file.");
}
int filter_id = std::stoi(get_node_value(node_filt, "id"));
// Check to make sure 'id' hasn't been used
if (model::filter_map.find(filter_id) != model::filter_map.end()) {
fatal_error("Two or more filters use the same unique ID: "
+ std::to_string(filter_id));
}
// Convert filter type to lower case
std::string s;
if (check_for_node(node_filt, "type")) {
s = get_node_value(node_filt, "type", true);
}
// Allocate according to the filter type
Filter* f = allocate_filter(s);
// Read filter data from XML
f->from_xml(node_filt);
// Set filter id
f->id_ = filter_id;
model::filter_map[filter_id] = model::tally_filters.size() - 1;
// Initialize filter
f->initialize();
auto f = Filter::create(node_filt);
}
// ==========================================================================
@ -657,229 +892,7 @@ void read_tallies_xml()
}
for (auto node_tal : root.children("tally")) {
model::tallies.push_back(std::make_unique<Tally>());
auto& t {model::tallies.back()};
t->init_from_xml(node_tal);
// Copy and set tally id
if (!check_for_node(node_tal, "id")) {
fatal_error("Must specify id for tally in tally XML file.");
}
t->id_ = std::stoi(get_node_value(node_tal, "id"));
model::tally_map[t->id_] = model::tallies.size() - 1;
// Copy tally name
if (check_for_node(node_tal, "name")) {
t->name_ = get_node_value(node_tal, "name");
}
// =======================================================================
// READ DATA FOR FILTERS
// Check if user is using old XML format and throw an error if so
if (check_for_node(node_tal, "filter")) {
fatal_error("Tally filters must be specified independently of "
"tallies in a <filter> element. The <tally> element itself should "
"have a list of filters that apply, e.g., <filters>1 2</filters> "
"where 1 and 2 are the IDs of filters specified outside of "
"<tally>.");
}
// Determine number of filters
std::vector<int> filters;
if (check_for_node(node_tal, "filters")) {
filters = get_node_array<int>(node_tal, "filters");
}
// Allocate and store filter user ids
if (!filters.empty()) {
std::vector<int> filter_indices;
for (int filter_id : filters) {
// Determine if filter ID is valid
auto it = model::filter_map.find(filter_id);
if (it == model::filter_map.end()) {
fatal_error("Could not find filter " + std::to_string(filter_id)
+ " specified on tally " + std::to_string(t->id_));
}
// Store the index of the filter
filter_indices.push_back(it->second);
}
// Set the filters
t->set_filters(filter_indices.data(), filter_indices.size());
}
// Check for the presence of certain filter types
bool has_energyout = t->energyout_filter_ >= 0;
int particle_filter_index = C_NONE;
for (int j = 0; j < t->filters().size(); ++j) {
int i_filter = t->filters(j);
const auto& f = model::tally_filters[i_filter].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
if (pf) particle_filter_index = i_filter;
// Change the tally estimator if a filter demands it
std::string filt_type = f->type();
if (filt_type == "energyout" || filt_type == "legendre") {
t->estimator_ = ESTIMATOR_ANALOG;
} else if (filt_type == "sphericalharmonics") {
auto sf = dynamic_cast<SphericalHarmonicsFilter*>(f);
if (sf->cosine_ == SphericalHarmonicsCosine::scatter) {
t->estimator_ = ESTIMATOR_ANALOG;
}
} else if (filt_type == "spatiallegendre" || filt_type == "zernike"
|| filt_type == "zernikeradial") {
t->estimator_ = ESTIMATOR_COLLISION;
}
}
// =======================================================================
// READ DATA FOR NUCLIDES
t->set_nuclides(node_tal);
// =======================================================================
// READ DATA FOR SCORES
t->set_scores(node_tal);
if (!check_for_node(node_tal, "scores")) {
fatal_error("No scores specified on tally " + std::to_string(t->id_)
+ ".");
}
// Check if tally is compatible with particle type
if (settings::photon_transport) {
if (particle_filter_index == C_NONE) {
for (int score : t->scores_) {
switch (score) {
case SCORE_INVERSE_VELOCITY:
fatal_error("Particle filter must be used with photon "
"transport on and inverse velocity score");
break;
case SCORE_FLUX:
case SCORE_TOTAL:
case SCORE_SCATTER:
case SCORE_NU_SCATTER:
case SCORE_ABSORPTION:
case SCORE_FISSION:
case SCORE_NU_FISSION:
case SCORE_CURRENT:
case SCORE_EVENTS:
case SCORE_DELAYED_NU_FISSION:
case SCORE_PROMPT_NU_FISSION:
case SCORE_DECAY_RATE:
warning("Particle filter is not used with photon transport"
" on and " + reaction_name(score) + " score.");
break;
}
}
} else {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (auto p : pf->particles_) {
if (p == Particle::Type::electron ||
p == Particle::Type::positron) {
t->estimator_ = ESTIMATOR_ANALOG;
}
}
}
} else {
if (particle_filter_index >= 0) {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (auto p : pf->particles_) {
if (p != Particle::Type::neutron) {
warning("Particle filter other than NEUTRON used with photon "
"transport turned off. All tallies for particle type " +
std::to_string(static_cast<int>(p)) + " will have no scores");
}
}
}
}
// Check for a tally derivative.
if (check_for_node(node_tal, "derivative")) {
int deriv_id = std::stoi(get_node_value(node_tal, "derivative"));
// Find the derivative with the given id, and store it's index.
auto it = model::tally_deriv_map.find(deriv_id);
if (it == model::tally_deriv_map.end()) {
fatal_error("Could not find derivative " + std::to_string(deriv_id)
+ " specified on tally " + std::to_string(t->id_));
}
t->deriv_ = it->second;
// Only analog or collision estimators are supported for differential
// tallies.
if (t->estimator_ == ESTIMATOR_TRACKLENGTH) {
t->estimator_ = ESTIMATOR_COLLISION;
}
const auto& deriv = model::tally_derivs[t->deriv_];
if (deriv.variable == DIFF_NUCLIDE_DENSITY
|| deriv.variable == DIFF_TEMPERATURE) {
for (int i_nuc : t->nuclides_) {
if (has_energyout && i_nuc == -1) {
fatal_error("Error on tally " + std::to_string(t->id_)
+ ": Cannot use a 'nuclide_density' or 'temperature' "
"derivative on a tally with an outgoing energy filter and "
"'total' nuclide rate. Instead, tally each nuclide in the "
"material individually.");
// Note that diff tallies with these characteristics would work
// correctly if no tally events occur in the perturbed material
// (e.g. pertrubing moderator but only tallying fuel), but this
// case would be hard to check for by only reading inputs.
}
}
}
}
// If settings.xml trigger is turned on, create tally triggers
if (settings::trigger_on) {
t->init_triggers(node_tal);
}
// =======================================================================
// SET TALLY ESTIMATOR
// Check if user specified estimator
if (check_for_node(node_tal, "estimator")) {
std::string est = get_node_value(node_tal, "estimator");
if (est == "analog") {
t->estimator_ = ESTIMATOR_ANALOG;
} else if (est == "tracklength" || est == "track-length"
|| est == "pathlength" || est == "path-length") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (t->estimator_ == ESTIMATOR_ANALOG) {
fatal_error("Cannot use track-length estimator for tally "
+ std::to_string(t->id_));
}
// Set estimator to track-length estimator
t->estimator_ = ESTIMATOR_TRACKLENGTH;
} else if (est == "collision") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (t->estimator_ == ESTIMATOR_ANALOG) {
fatal_error("Cannot use collision estimator for tally " +
std::to_string(t->id_));
}
// Set estimator to collision estimator
t->estimator_ = ESTIMATOR_COLLISION;
} else {
fatal_error("Invalid estimator '" + est + "' on tally " +
std::to_string(t->id_));
}
}
model::tallies.push_back(std::make_unique<Tally>(node_tal));
}
}
@ -1057,7 +1070,7 @@ openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end)
if (index_start) *index_start = model::tallies.size();
if (index_end) *index_end = model::tallies.size() + n - 1;
for (int i = 0; i < n; ++i) {
model::tallies.push_back(std::make_unique<Tally>());
model::tallies.push_back(std::make_unique<Tally>(-1));
}
return 0;
}
@ -1141,14 +1154,7 @@ openmc_tally_set_id(int32_t index, int32_t id)
return OPENMC_E_OUT_OF_BOUNDS;
}
if (model::tally_map.find(id) != model::tally_map.end()) {
set_errmsg("Two or more tallies use the same unique ID: "
+ std::to_string(id));
return OPENMC_E_INVALID_ID;
}
model::tallies[index]->id_ = id;
model::tally_map[id] = index;
model::tallies[index]->set_id(id);
return 0;
}
@ -1288,7 +1294,7 @@ openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides)
}
extern "C" int
openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n)
openmc_tally_get_filters(int32_t index, const int32_t** indices, size_t* n)
{
if (index < 0 || index >= model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
@ -1301,7 +1307,7 @@ openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n)
}
extern "C" int
openmc_tally_set_filters(int32_t index, int n, const int32_t* indices)
openmc_tally_set_filters(int32_t index, size_t n, const int32_t* indices)
{
// Make sure the index fits in the array bounds.
if (index < 0 || index >= model::tallies.size()) {
@ -1311,9 +1317,15 @@ openmc_tally_set_filters(int32_t index, int n, const int32_t* indices)
// Set the filters.
try {
model::tallies[index]->set_filters(indices, n);
// Convert indices to filter pointers
std::vector<Filter*> filters;
for (gsl::index i = 0; i < n; ++i) {
int32_t i_filt = indices[i];
filters.push_back(model::tally_filters.at(i_filt).get());
}
model::tallies[index]->set_filters(filters);
} catch (const std::out_of_range& ex) {
set_errmsg(ex.what());
set_errmsg("Index in tally filter array out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}

72
src/tallies/tally_scoring.cpp
View file

@ -69,7 +69,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end)
if (!match.bins_present_) {
match.bins_.clear();
match.weights_.clear();
for (auto i = 0; i < model::tally_filters[i_filt]->n_bins_; ++i) {
for (auto i = 0; i < model::tally_filters[i_filt]->n_bins(); ++i) {
match.bins_.push_back(i);
match.weights_.push_back(1.0);
}
@ -209,14 +209,14 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
if (tally.deriv_ != C_NONE)
apply_derivative_to_score(p, i_tally, 0, 0., SCORE_NU_FISSION, score);
if (!settings::run_CE && eo_filt.matches_transport_groups_) {
if (!settings::run_CE && eo_filt.matches_transport_groups()) {
// determine outgoing energy group from fission bank
auto g_out = static_cast<int>(bank.E);
// modify the value so that g_out = 1 corresponds to the highest energy
// bin
g_out = eo_filt.n_bins_ - g_out;
g_out = eo_filt.n_bins() - g_out;
// change outgoing energy bin
simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out;
@ -231,11 +231,11 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
}
// Set EnergyoutFilter bin index
if (E_out < eo_filt.bins_.front() || E_out > eo_filt.bins_.back()) {
if (E_out < eo_filt.bins().front() || E_out > eo_filt.bins().back()) {
continue;
} else {
auto i_match = lower_bound_index(eo_filt.bins_.begin(),
eo_filt.bins_.end(), E_out);
auto i_match = lower_bound_index(eo_filt.bins().begin(),
eo_filt.bins().end(), E_out);
simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match;
}
@ -272,8 +272,8 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
model::tally_filters[i_dg_filt].get())};
// Loop over delayed group bins until the corresponding bin is found
for (auto d_bin = 0; d_bin < dg_filt.n_bins_; ++d_bin) {
if (dg_filt.groups_[d_bin] == g) {
for (auto d_bin = 0; d_bin < dg_filt.n_bins(); ++d_bin) {
if (dg_filt.groups()[d_bin] == g) {
// Find the filter index and weight for this filter combination
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
@ -632,8 +632,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield = data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = p->wgt_absorb_ * yield
@ -670,8 +670,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
score = simulation::keff * p->wgt_bank_ / p->n_bank_
* p->n_delayed_bank_[d-1] * flux;
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
@ -693,8 +693,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield = data::nuclides[i_nuclide]
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = p->neutron_xs_[i_nuclide].fission * yield
@ -722,8 +722,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield = data::nuclides[j_nuclide]
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = p->neutron_xs_[j_nuclide].fission * yield
@ -770,8 +770,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield
= nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = rxn.products_[d].decay_rate_;
@ -830,8 +830,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Find the corresponding filter bin and then score
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
if (d == g)
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
@ -852,8 +852,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield
= nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = rxn.products_[d].decay_rate_;
@ -892,8 +892,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
if (nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
auto yield
= nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = rxn.products_[d].decay_rate_;
@ -1757,8 +1757,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *=
@ -1807,8 +1807,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
score = simulation::keff * p->wgt_bank_ / p->n_bank_
* p->n_delayed_bank_[d-1] * flux;
if (i_nuclide >= 0) {
@ -1839,8 +1839,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
if (i_nuclide >= 0) {
score = flux * atom_density
* get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION,
@ -1879,8 +1879,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *=
@ -1959,8 +1959,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Find the corresponding filter bin and then score
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
if (d == g)
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
@ -1978,8 +1978,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin];
if (i_nuclide >= 0) {
score += atom_density * flux
* get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,

3
tests/regression_tests/complex_cell/test.py
View file

@ -1,5 +1,8 @@
from tests.testing_harness import TestHarness
import sys
import openmc.capi
def test_complex_cell():
harness = TestHarness('statepoint.10.h5')

BIN
tests/regression_tests/dagmc/legacy/dagmc.h5m
View file

Binary file not shown.

10
tests/regression_tests/dagmc/legacy/test.py
View file

@ -1,6 +1,5 @@
import openmc
import openmc.capi
from openmc.stats import Box
import pytest
from tests.testing_harness import PyAPITestHarness
@ -17,8 +16,10 @@ def test_dagmc():
model.settings.inactive = 0
model.settings.particles = 100
source = openmc.Source(space=Box([-4, -4, -4],
[ 4, 4, 4]))
source_box = openmc.stats.Box([-4, -4, -4],
[ 4, 4, 4])
source = openmc.Source(space=source_box)
model.settings.source = source
model.settings.dagmc = True
@ -45,5 +46,4 @@ def test_dagmc():
mats = openmc.Materials([u235, water])
model.materials = mats
harness = PyAPITestHarness('statepoint.5.h5', model=model)
harness.main()
model.export_to_xml()

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