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

Merge remote-tracking branch 'mitcrpg/develop' into multipole-final

Browse source
This commit is contained in:
Colin Josey 2016-02-04 20:30:20 -05:00
commit b9504fd542
226 changed files with 3470 additions and 8609 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.gitignore vendored
View file

@ -43,6 +43,8 @@ results_test.dat
# Test build files
tests/build/
tests/coverage/
tests/memcheck/
tests/ctestscript.run
# HDF5 files
@ -60,6 +62,7 @@ data/nndc
#Images
*.ppm
*.voxel
# PyCharm project configuration files
.idea

109
CMakeLists.txt
View file

@ -126,8 +126,8 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(APPEND f90flags -O3)
list(APPEND cflags -O3)
list(APPEND f90flags -O3 -flto -fuse-linker-plugin)
list(APPEND cflags -O3 -flto -fuse-linker-plugin)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
@ -335,113 +335,36 @@ include(CTest)
# Get a list of all the tests to run
file(GLOB_RECURSE TESTS ${CMAKE_CURRENT_SOURCE_DIR}/tests/test_*.py)
# Check for MEM_CHECK and COVERAGE variables
if (DEFINED ENV{MEM_CHECK})
set(MEM_CHECK $ENV{MEM_CHECK})
else(DEFINED ENV{MEM_CHECK})
set(MEM_CHECK FALSE)
endif(DEFINED ENV{MEM_CHECK})
if (DEFINED ENV{COVERAGE})
set(COVERAGE $ENV{COVERAGE})
else(DEFINED ENV{COVERAGE})
set(COVERAGE FALSE)
endif(DEFINED ENV{COVERAGE})
# Loop through all the tests
foreach(test ${TESTS})
# Get test information
get_filename_component(TEST_NAME ${test} NAME)
get_filename_component(TEST_PATH ${test} PATH)
# Check for running standard tests (no valgrind, no gcov)
if(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
if (DEFINED ENV{MEM_CHECK})
# Generate input files if needed
if (NOT EXISTS "${TEST_PATH}/geometry.xml")
execute_process(COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --build-inputs
WORKING_DIRECTORY ${TEST_PATH})
endif()
# Add serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc>)
else()
# Check serial/parallel
if (${MPI_ENABLED})
# Preform a parallel test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>
--mpi_exec $ENV{MPI_DIR}/bin/mpiexec)
else(${MPI_ENABLED})
else()
# Perform a serial test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND ${PYTHON_EXECUTABLE} ${TEST_NAME} --exe $<TARGET_FILE:openmc>)
endif(${MPI_ENABLED})
# Handle special case for valgrind and gcov (run openmc directly, no python)
else(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
# If a plot test is encountered, run with "-p"
if (${test} MATCHES "test_plot")
# Perform serial valgrind and coverage test with plot flag
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> -p ${TEST_PATH})
elseif(${test} MATCHES "test_filter_distribcell")
# Add each case for distribcell tests
add_test(NAME ${TEST_NAME}_case-1
WORKING_DIRECTORY ${TEST_PATH}/case-1
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-1)
add_test(NAME ${TEST_NAME}_case-2
WORKING_DIRECTORY ${TEST_PATH}/case-2
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-2)
add_test(NAME ${TEST_NAME}_case-3
WORKING_DIRECTORY ${TEST_PATH}/case-3
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-3)
add_test(NAME ${TEST_NAME}_case-4
WORKING_DIRECTORY ${TEST_PATH}/case-4
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH}/case-4)
# If a restart test is encounted, need to run with -r and restart file(s)
elseif(${test} MATCHES "restart")
# Handle restart tests separately
if(${test} MATCHES "test_statepoint_restart")
set(RESTART_FILE statepoint.07.h5)
elseif(${test} MATCHES "test_sourcepoint_restart")
set(RESTART_FILE statepoint.07.h5 source.07.h5)
elseif(${test} MATCHES "test_particle_restart_eigval")
set(RESTART_FILE particle_9_555.h5)
elseif(${test} MATCHES "test_particle_restart_fixed")
set(RESTART_FILE particle_7_928.h5)
else(${test} MATCHES "test_statepoint_restart")
message(FATAL_ERROR "Restart test ${test} not recognized")
endif(${test} MATCHES "test_statepoint_restart")
# Perform serial valgrind and coverage test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH})
# Perform serial valgrind and coverage restart test
add_test(NAME ${TEST_NAME}_restart
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> -r ${RESTART_FILE} ${TEST_PATH})
# Set test dependency
set_tests_properties(${TEST_NAME}_restart PROPERTIES DEPENDS ${TEST_NAME})
# Handle standard tests for valgrind and gcov
else(${test} MATCHES "test_plot")
# Perform serial valgrind and coverage test
add_test(NAME ${TEST_NAME}
WORKING_DIRECTORY ${TEST_PATH}
COMMAND $<TARGET_FILE:openmc> ${TEST_PATH})
endif(${test} MATCHES "test_plot")
endif(NOT ${MEM_CHECK} AND NOT ${COVERAGE})
endif()
endif()
endforeach(test)

10
docs/source/_static/theme_overrides.css Normal file
View file

@ -0,0 +1,10 @@
/* override table width restrictions */
.wy-table-responsive table td, .wy-table-responsive table th {
white-space: normal;
}
.wy-table-responsive {
margin-bottom: 24px;
max-width: 100%;
overflow: visible;
}

5
docs/source/conf.py
View file

@ -152,7 +152,10 @@ html_title = "OpenMC Documentation"
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
#html_static_path = ['_static']
html_static_path = ['_static']
def setup(app):
app.add_stylesheet('theme_overrides.css')
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.

2
docs/source/publications.rst
View file

@ -28,7 +28,7 @@ Benchmarking
- Khurrum S. Chaudri and Sikander M. Mirza, "Burnup dependent Monte Carlo
neutron physics calculations of IAEA MTR benchmark," *Prog. Nucl. Energy*,
**81**, 43-52 (2015). `<http://dx.doi.org/j.pnucene.2014.12.018>`_
**81**, 43-52 (2015). `<http://dx.doi.org/10.1016/j.pnucene.2014.12.018>`_
- Daniel J. Kelly, Brian N. Aviles, Paul K. Romano, Bryan R. Herman,
Nicholas E. Horelik, and Benoit Forget, "Analysis of select BEAVRS PWR

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

File diff suppressed because one or more lines are too long

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

@ -1452,7 +1452,6 @@
],
"source": [
"# Generate tracks for OpenMOC\n",
"openmoc_geometry.initializeFlatSourceRegions()\n",
"track_generator = openmoc.TrackGenerator(openmoc_geometry, num_azim=32, spacing=0.1)\n",
"track_generator.generateTracks()\n",
"\n",
@ -1638,7 +1637,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10"
"version": "2.7.6"
}
},
"nbformat": 4,

44
docs/source/quickinstall.rst
View file

@ -12,24 +12,55 @@ OpenMC, see :ref:`usersguide_install` in the User's Manual.
Installing on Ubuntu through PPA
--------------------------------
For users with Ubuntu 11.10 or later, a binary package for OpenMC is available
through a `Personal Package Archive`_ (PPA) and can be installed through the `APT
package manager`_. Simply enter the following commands into the terminal:
For users with Ubuntu 15.04 or later, a binary package for OpenMC is available
through a `Personal Package Archive`_ (PPA) and can be installed through the
`APT package manager`_. First, add the following PPA to the repository sources:
.. code-block:: sh
sudo apt-add-repository ppa:paulromano/staging
Next, resynchronize the package index files:
.. code-block:: sh
sudo apt-get update
Now OpenMC should be recognized within the repository and can be installed:
.. code-block:: sh
sudo apt-get install openmc
Currently, the binary package does not allow for parallel simulations or use of
HDF5_. Users who need such capabilities should build OpenMC from source as is
described in :ref:`usersguide_install`.
Binary packages from this PPA may exist for earlier versions of Ubuntu, but they
are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
.. _HDF5: http://www.hdfgroup.org/HDF5/
---------------------------------------
Installing from Source on Ubuntu 15.04+
---------------------------------------
To build OpenMC from source, several :ref:`prerequisites <prerequisites>` are
needed. If you are using Ubuntu 15.04 or higher, all prerequisites can be
installed directly from the package manager.
.. code-block:: sh
sudo apt-get install gfortran
sudo apt-get install cmake
sudo apt-get install libhdf5-dev
After the packages have been installed, follow the instructions below for
building and installing OpenMC from source.
.. note:: Before Ubuntu 15.04, the HDF5 package included in the Ubuntu Package
archive was not built with support for the Fortran 2003 HDF5
interface, which is needed by OpenMC. If you are using Ubuntu 14.10 or
before you will need to build HDF5 from source.
-------------------------------------------
Installing from Source on Linux or Mac OS X
-------------------------------------------
@ -42,7 +73,6 @@ entering the following commands in a terminal:
git clone https://github.com/mit-crpg/openmc.git
cd openmc
git checkout -b master origin/master
mkdir build && cd build
cmake ..
make

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

@ -1484,114 +1484,203 @@ The ``<tally>`` element accepts the following sub-elements:
*Default*: ``tracklength`` but will revert to ``analog`` if necessary.
:scores:
A space-separated list of the desired responses to be accumulated. Accepted
options are "flux", "total", "scatter", "absorption", "fission",
"nu-fission", "delayed-nu-fission", "kappa-fission", "nu-scatter",
"scatter-N", "scatter-PN", "scatter-YN", "nu-scatter-N", "nu-scatter-PN",
"nu-scatter-YN", "flux-YN", "total-YN", "current", "inverse-velocity" and
"events". These correspond to the following physical quantities:
A space-separated list of the desired responses to be accumulated. The accepted
options are listed in the following tables:
:flux:
Total flux in particle-cm per source particle.
.. table:: **Flux scores: units are particle-cm per source particle.**
.. note::
The ``analog`` estimator is actually identical to the ``collision``
estimator for the flux score.
+----------------------+---------------------------------------------------+
|Score | Description |
+======================+===================================================+
|flux |Total flux. |
+----------------------+---------------------------------------------------+
|flux-YN |Spherical harmonic expansion of the direction of |
| |motion :math:`\left(\Omega\right)` of the total |
| |flux. This score will tally all of the harmonic |
| |moments of order 0 to N. N must be between 0 and |
| |10. |
+----------------------+---------------------------------------------------+
:total:
Total reaction rate in reactions per source particle.
.. table:: **Reaction scores: units are reactions per source particle.**
:scatter:
Total scattering rate. Can also be identified with the ``scatter-0``
response type. Units are reactions per source particle.
+----------------------+---------------------------------------------------+
|Score | Description |
+======================+===================================================+
|absorption |Total absorption rate. This accounts for all |
| |reactions which do not produce secondary neutrons. |
+----------------------+---------------------------------------------------+
|elastic |Elastic scattering reaction rate. |
+----------------------+---------------------------------------------------+
|fission |Total fission reaction rate. |
+----------------------+---------------------------------------------------+
|scatter |Total scattering rate. Can also be identified with |
| |the "scatter-0" response type. |
+----------------------+---------------------------------------------------+
|scatter-N |Tally the N\ :sup:`th` \ scattering moment, where N|
| |is the Legendre expansion order of the change in |
| |particle angle :math:`\left(\mu\right)`. N must be |
| |between 0 and 10. As an example, tallying the 2\ |
| |:sup:`nd` \ scattering moment would be specified as|
| |``<scores>scatter-2</scores>``. |
+----------------------+---------------------------------------------------+
|scatter-PN |Tally all of the scattering moments from order 0 to|
| |N, where N is the Legendre expansion order of the |
| |change in particle angle |
| |:math:`\left(\mu\right)`. That is, "scatter-P1" is |
| |equivalent to requesting tallies of "scatter-0" and|
| |"scatter-1". Like for "scatter-N", N must be |
| |between 0 and 10. As an example, tallying up to the|
| |2\ :sup:`nd` \ scattering moment would be specified|
| |as ``<scores> scatter-P2 </scores>``. |
+----------------------+---------------------------------------------------+
|scatter-YN |"scatter-YN" is similar to "scatter-PN" except an |
| |additional expansion is performed for the incoming |
| |particle direction :math:`\left(\Omega\right)` |
| |using the real spherical harmonics. This is useful|
| |for performing angular flux moment weighting of the|
| |scattering moments. Like "scatter-PN", "scatter-YN"|
| |will tally all of the moments from order 0 to N; N |
| |again must be between 0 and 10. |
+----------------------+---------------------------------------------------+
|total |Total reaction rate. |
+----------------------+---------------------------------------------------+
|total-YN |The total reaction rate expanded via spherical |
| |harmonics about the direction of motion of the |
| |neutron, :math:`\Omega`. This score will tally all |
| |of the harmonic moments of order 0 to N. N must be|
| |between 0 and 10. |
+----------------------+---------------------------------------------------+
|(n,2nd) |(n,2nd) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2n) |(n,2n) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,3n) |(n,3n) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,na) |(n,n\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,n3a) |(n,n3\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2na) |(n,2n\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,3na) |(n,3n\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,np) |(n,np) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,n2a) |(n,n2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2n2a) |(n,2n2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,nd) |(n,nd) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,nt) |(n,nt) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,nHe-3) |(n,n\ :sup:`3`\ He) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,nd2a) |(n,nd2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,nt2a) |(n,nt2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,4n) |(n,4n) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2np) |(n,2np) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,3np) |(n,3np) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,n2p) |(n,n2p) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,n*X*) |Level inelastic scattering reaction rate. The *X* |
| |indicates what which inelastic level, e.g., (n,n3) |
| |is third-level inelastic scattering. |
+----------------------+---------------------------------------------------+
|(n,nc) |Continuum level inelastic scattering reaction rate.|
+----------------------+---------------------------------------------------+
|(n,gamma) |Radiative capture reaction rate. |
+----------------------+---------------------------------------------------+
|(n,p) |(n,p) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,d) |(n,d) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,t) |(n,t) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,3He) |(n,\ :sup:`3`\ He) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,a) |(n,\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2a) |(n,2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,3a) |(n,3\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2p) |(n,2p) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,pa) |(n,p\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,t2a) |(n,t2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,d2a) |(n,d2\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,pd) |(n,pd) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,pt) |(n,pt) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,da) |(n,d\ :math:`\alpha`\ ) reaction rate. |
+----------------------+---------------------------------------------------+
|*Arbitrary integer* |An arbitrary integer is interpreted to mean the |
| |reaction rate for a reaction with a given ENDF MT |
| |number. |
+----------------------+---------------------------------------------------+
:absorption:
Total absorption rate. This accounts for all reactions which do not
produce secondary neutrons. Units are reactions per source particle.
.. table:: **Particle production scores: units are particles produced per
source particles.**
:fission:
Total fission rate in reactions per source particle.
+----------------------+---------------------------------------------------+
|Score | Description |
+======================+===================================================+
|delayed-nu-fission |Total production of delayed neutrons due to |
| |fission. |
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
|nu-scatter, |These scores are similar in functionality to their |
|nu-scatter-N, |``scatter*`` equivalents except the total |
|nu-scatter-PN, |production of neutrons due to scattering is scored |
|nu-scatter-YN |vice simply the scattering rate. This accounts for |
| |multiplicity from (n,2n), (n,3n), and (n,4n) |
| |reactions. |
+----------------------+---------------------------------------------------+
:nu-fission:
Total production of neutrons due to fission. Units are neutrons produced
per source neutron.
.. table:: **Miscellaneous scores: units are indicated for each.**
:delayed-nu-fission:
Total production of delayed neutrons due to fission. Units are neutrons produced
per source neutron.
+----------------------+---------------------------------------------------+
|Score | Description |
+======================+===================================================+
|current |Partial currents on the boundaries of each cell in |
| |a mesh. Units are particles per source |
| |particle. Note that this score can only be used if |
| |a mesh filter has been specified. Furthermore, it |
| |may not be used in conjunction with any other |
| |score. |
+----------------------+---------------------------------------------------+
|events |Number of scoring events. Units are events per |
| |source particle. |
+----------------------+---------------------------------------------------+
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
| |fission product kinetic energy, prompt and delayed |
| |neutron kinetic energies, prompt and delayed |
| |:math:`\gamma`-ray total energies, and the total |
| |energy released by the delayed :math:`\beta` |
| |particles. The neutrino energy does not contribute |
| |to this response. The prompt and delayed |
| |:math:`\gamma`-rays are assumed to deposit their |
| |energy locally. Units are MeV per source particle. |
+----------------------+---------------------------------------------------+
:kappa-fission:
The recoverable energy production rate due to fission. The recoverable
energy is defined as the fission product kinetic energy, prompt and
delayed neutron kinetic energies, prompt and delayed :math:`\gamma`-ray
total energies, and the total energy released by the delayed :math:`\beta`
particles. The neutrino energy does not contribute to this response. The
prompt and delayed :math:`\gamma`-rays are assumed to deposit their energy
locally. Units are MeV per source particle.
:scatter-N:
Tally the N\ :sup:`th` \ scattering moment, where N is the Legendre
expansion order of the change in particle angle :math:`\left(\mu\right)`.
N must be between 0 and 10. As an example, tallying the 2\ :sup:`nd` \
scattering moment would be specified as ``<scores> scatter-2
</scores>``. Units are reactions per source particle.
:scatter-PN:
Tally all of the scattering moments from order 0 to N, where N is the
Legendre expansion order of the change in particle angle
:math:`\left(\mu\right)`. That is, ``scatter-P1`` is equivalent to
requesting tallies of ``scatter-0`` and ``scatter-1``. Like for
``scatter-N``, N must be between 0 and 10. As an example, tallying up to
the 2\ :sup:`nd` \ scattering moment would be specified as ``<scores>
scatter-P2 </scores>``. Units are reactions per source particle.
:scatter-YN:
``scatter-YN`` is similar to ``scatter-PN`` except an additional expansion
is performed for the incoming particle direction
:math:`\left(\Omega\right)` using the real spherical harmonics. This is
useful for performing angular flux moment weighting of the scattering
moments. Like ``scatter-PN``, ``scatter-YN`` will tally all of the moments
from order 0 to N; N again must be between 0 and 10. Units are reactions
per source particle.
:nu-scatter, nu-scatter-N, nu-scatter-PN, nu-scatter-YN:
These scores are similar in functionality to their ``scatter*``
equivalents except the total production of neutrons due to scattering is
scored vice simply the scattering rate. This accounts for multiplicity
from (n,2n), (n,3n), and (n,4n) reactions. Units are neutrons produced per
source particle.
:flux-YN:
Spherical harmonic expansion of the direction of motion
:math:`\left(\Omega\right)` of the total flux. This score will tally all
of the harmonic moments of order 0 to N. N must be between 0
and 10. Units are particle-cm per source particle.
:total-YN:
The total reaction rate expanded via spherical harmonics about the
direction of motion of the neutron, :math:`\Omega`.
This score will tally all of the harmonic moments of order 0 to N. N must
be between 0 and 10. Units are reactions per source particle.
:current:
Partial currents on the boundaries of each cell in a mesh. Units are
particles per source particle.
.. note::
This score can only be used if a mesh filter has been
specified. Furthermore, it may not be used in conjunction with any
other score.
:inverse-velocity:
The flux-weighted inverse velocity where the velocity is in units of
centimeters per second.
.. note::
The ``analog`` estimator is actually identical to the ``collision``
estimator for the inverse-velocity score.
:events:
Number of scoring events. Units are events per source particle.
.. note::
The ``analog`` estimator is actually identical to the ``collision``
estimator for the flux and inverse-velocity scores.
:trigger:
Precision trigger applied to all filter bins and nuclides for this tally.

9
docs/source/usersguide/install.rst
View file

@ -9,8 +9,8 @@ Installing on Ubuntu with PPA
-----------------------------
For users with Ubuntu 15.04 or later, a binary package for OpenMC is available
through a Personal Package Archive (PPA) and can be installed through the APT
package manager. First, add the following PPA to the repository sources:
through a `Personal Package Archive`_ (PPA) and can be installed through the
`APT package manager`_. First, add the following PPA to the repository sources:
.. code-block:: sh
@ -31,10 +31,15 @@ Now OpenMC should be recognized within the repository and can be installed:
Binary packages from this PPA may exist for earlier versions of Ubuntu, but they
are no longer supported.
.. _Personal Package Archive: https://launchpad.net/~paulromano/+archive/staging
.. _APT package manager: https://help.ubuntu.com/community/AptGet/Howto
--------------------
Building from Source
--------------------
.. _prerequisites:
Prerequisites
-------------

409
openmc/aggregate.py
View file

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

423
openmc/cross.py → openmc/arithmetic.py
View file

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

4
openmc/filter.py
View file

@ -673,9 +673,11 @@ class Filter(object):
# Assign entry to Lattice Multi-index column
else:
# Reverse y index per lattice ordering in OpenCG
level_dict[lat_id_key][offset] = coords._lattice._id
level_dict[lat_x_key][offset] = coords._lat_x
level_dict[lat_y_key][offset] = coords._lat_y
level_dict[lat_y_key][offset] = \
coords._lattice.dimension[1] - coords._lat_y - 1
level_dict[lat_z_key][offset] = coords._lat_z
# Move to next node in LocalCoords linked list

7
openmc/mgxs/library.py
View file

@ -396,7 +396,7 @@ class Library(object):
cv.check_type('statepoint', statepoint, openmc.StatePoint)
if not statepoint.with_summary:
if statepoint.summary is None:
msg = 'Unable to load data from a statepoint which has not been ' \
'linked with a summary file'
raise ValueError(msg)
@ -568,8 +568,9 @@ class Library(object):
for domain in self.domains:
for mgxs_type in self.mgxs_types:
mgxs = subdomain_avg_library.get_mgxs(domain, mgxs_type)
avg_mgxs = mgxs.get_subdomain_avg_xs()
subdomain_avg_library.all_mgxs[domain.id][mgxs_type] = avg_mgxs
if mgxs.domain_type == 'distribcell':
avg_mgxs = mgxs.get_subdomain_avg_xs()
subdomain_avg_library.all_mgxs[domain.id][mgxs_type] = avg_mgxs
return subdomain_avg_library

81
openmc/mgxs/mgxs.py
View file

@ -583,7 +583,7 @@ class MGXS(object):
cv.check_type('statepoint', statepoint, openmc.statepoint.StatePoint)
if not statepoint.with_summary:
if statepoint.summary is None:
msg = 'Unable to load data from a statepoint which has not been ' \
'linked with a summary file'
raise ValueError(msg)
@ -612,6 +612,11 @@ class MGXS(object):
filters = []
filter_bins = []
# Clear any tallies previously loaded from a statepoint
self._tallies = None
self._xs_tally = None
self._rxn_rate_tally = None
# Find, slice and store Tallies from StatePoint
# The tally slicing is needed if tally merging was used
for tally_type, tally in self.tallies.items():
@ -673,14 +678,14 @@ class MGXS(object):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append('energy')
@ -838,55 +843,25 @@ class MGXS(object):
"""
# Construct a collection of the subdomain filter bins to average across
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains)
else:
subdomains = [0]
subdomains = None
# Clone this MGXS to initialize the subdomain-averaged version
avg_xs = copy.deepcopy(self)
avg_xs._rxn_rate_tally = None
avg_xs._xs_tally = None
avg_xs._sparse = False
# If domain is distribcell, make the new domain 'cell'
if self.domain_type == 'distribcell':
avg_xs.domain_type = 'cell'
# Average each of the tallies across subdomains
for tally_type, tally in avg_xs.tallies.items():
tally_avg = tally.summation(filter_type=self.domain_type,
filter_bins=subdomains)
avg_xs.tallies[tally_type] = tally_avg
# Make condensed tally derived and null out sum, sum_sq
tally._derived = True
tally._sum = None
tally._sum_sq = None
# Get tally data arrays reshaped with one dimension per filter
mean = tally.get_reshaped_data(value='mean')
std_dev = tally.get_reshaped_data(value='std_dev')
# Get the mean, std. dev. across requested subdomains
mean = np.sum(mean[subdomains, ...], axis=0)
std_dev = np.sum(std_dev[subdomains, ...]**2, axis=0)
std_dev = np.sqrt(std_dev)
# If domain is distribcell, make subdomain-averaged a 'cell' domain
domain_filter = tally.find_filter(self._domain_type)
if domain_filter.type == 'distribcell':
domain_filter.type = 'cell'
domain_filter.num_bins = 1
# Reshape averaged data arrays with one dimension for all filters
mean = np.reshape(mean, tally.shape)
std_dev = np.reshape(std_dev, tally.shape)
# Override tally's data with the new condensed data
tally._mean = mean
tally._std_dev = std_dev
# Compute the subdomain-averaged multi-group cross section
avg_xs._domain_type = 'sum({0})'.format(self.domain_type)
avg_xs.sparse = self.sparse
return avg_xs
@ -911,7 +886,7 @@ class MGXS(object):
"""
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1040,7 +1015,7 @@ class MGXS(object):
xs_results = h5py.File(filename, 'w')
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1222,7 +1197,7 @@ class MGXS(object):
"""
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
if nuclides != 'all' and nuclides != 'sum':
cv.check_iterable_type('nuclides', nuclides, basestring)
@ -1249,7 +1224,7 @@ class MGXS(object):
df = self.xs_tally.get_pandas_dataframe(summary=summary)
# Remove the score column since it is homogeneous and redundant
if summary and self.domain_type == 'distribcell':
if summary and 'distribcell' in self.domain_type:
df = df.drop('score', level=0, axis=1)
else:
df = df.drop('score', axis=1)
@ -1282,7 +1257,7 @@ class MGXS(object):
columns = ['group in']
# Select out those groups the user requested
if groups != 'all':
if not isinstance(groups, basestring):
if 'group in' in df:
df = df[df['group in'].isin(groups)]
if 'group out' in df:
@ -1819,21 +1794,21 @@ class ScatterMatrixXS(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if in_groups != 'all':
if not isinstance(in_groups, basestring):
cv.check_iterable_type('groups', in_groups, Integral)
for group in in_groups:
filters.append('energy')
filter_bins.append((self.energy_groups.get_group_bounds(group),))
# Construct list of energy group bounds tuples for all requested groups
if out_groups != 'all':
if not isinstance(out_groups, basestring):
cv.check_iterable_type('groups', out_groups, Integral)
for group in out_groups:
filters.append('energyout')
@ -1917,7 +1892,7 @@ class ScatterMatrixXS(MGXS):
"""
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1956,12 +1931,6 @@ class ScatterMatrixXS(MGXS):
bounds = self.energy_groups.get_group_bounds(group)
string += template.format('', group, bounds[0], bounds[1])
if subdomains == 'all':
if self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
else:
subdomains = [self.domain.id]
# Loop over all subdomains
for subdomain in subdomains:
@ -2165,14 +2134,14 @@ class Chi(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append('energyout')

4
openmc/statepoint.py
View file

@ -449,10 +449,6 @@ class StatePoint(object):
def summary(self):
return self._summary
@property
def with_summary(self):
return False if self.summary is None else True
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)

4
openmc/summary.py
View file

@ -370,9 +370,9 @@ class Summary(object):
# Set the universes for the lattice
lattice.universes = universes
# Set the distribcell offsets for the lattice
if offsets is not None:
offsets = np.swapaxes(offsets, 0, 2)
lattice.offsets = offsets
lattice.offsets = offsets[:, ::-1, :]
# Add the Lattice to the global dictionary of all Lattices
self.lattices[index] = lattice

77
openmc/tallies.py
View file

@ -12,8 +12,7 @@ import sys
import numpy as np
from openmc import Mesh, Filter, Trigger, Nuclide
from openmc.cross import CrossScore, CrossNuclide, CrossFilter
from openmc.aggregate import AggregateScore, AggregateNuclide, AggregateFilter
from openmc.arithmetic import *
from openmc.filter import _FILTER_TYPES
import openmc.checkvalue as cv
from openmc.clean_xml import *
@ -1003,16 +1002,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins for the corresponding filter type in the filters
parameter. Each bins is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
for the cell instance ID for 'distribcell' Filters. Each bin is a
2-tuple of floats for 'energy' and 'energyout' filters corresponding
to the energy boundaries of the bin of interest. The bin is an
(x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell
of interest. The order of the bins in the list must correspond to
the filter_types parameter.
Returns
-------
@ -1174,16 +1174,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins for the corresponding filter type in the filters
parameter. Each bins is the integer ID for 'material', 'surface',
'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
for the cell instance ID for 'distribcell' Filters. Each bin is a
2-tuple of floats for 'energy' and 'energyout' filters corresponding
to the energy boundaries of the bin of interest. The bin is an
(x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell
of interest. The order of the bins in the list must correspond to
the filter_types parameter.
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
@ -2641,16 +2642,17 @@ class Tally(object):
A list of filter type strings
(e.g., ['mesh', 'energy']; default is [])
filter_bins : list of Iterables
A list of the filter bins corresponding to the filter_types
parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
in the list is the integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. Each bin is an integer for the
cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest. The order of the bins in the list must correspond to the
filter_types parameter.
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins to slice for the corresponding filter type in
the filters parameter. Each bins is the integer ID for 'material',
'surface', 'cell', 'cellborn', and 'universe' Filters. Each bin is
an integer for the cell instance ID for 'distribcell' Filters. Each
bin is a 2-tuple of floats for 'energy' and 'energyout' filters
corresponding to the energy boundaries of the bin of interest. The
bin is an (x,y,z) 3-tuple for 'mesh' filters corresponding to the
mesh cell of interest. The order of the bins in the list must
correspond to the filter_types parameter.
nuclides : list of str
A list of nuclide name strings
(e.g., ['U-235', 'U-238']; default is [])
@ -2677,11 +2679,11 @@ class Tally(object):
new_tally = copy.deepcopy(self)
new_tally.sparse = False
if self.sum is not None:
if not self.derived and self.sum is not None:
new_sum = self.get_values(scores, filters, filter_bins,
nuclides, 'sum')
new_tally.sum = new_sum
if self.sum_sq is not None:
if not self.derived and self.sum_sq is not None:
new_sum_sq = self.get_values(scores, filters, filter_bins,
nuclides, 'sum_sq')
new_tally.sum_sq = new_sum_sq
@ -2736,6 +2738,7 @@ class Tally(object):
for filter_bin in filter_bins[i]:
bin_index = find_filter.get_bin_index(filter_bin)
if filter_type in ['energy', 'energyout']:
bin_indices.extend([bin_index])
bin_indices.extend([bin_index, bin_index+1])
num_bins += 1
elif filter_type == 'distribcell':
@ -2745,7 +2748,7 @@ class Tally(object):
bin_indices.append(bin_index)
num_bins += 1
find_filter.bins = find_filter.bins[bin_indices]
find_filter.bins = np.unique(find_filter.bins[bin_indices])
find_filter.num_bins = num_bins
# Update the new tally's filter strides
@ -2951,10 +2954,10 @@ class Tally(object):
diag_indices[start:end] = indices + (i * new_filter.num_bins**2)
# Inject this Tally's data along the diagonal of the diagonalized Tally
if self.sum is not None:
if not self.derived and self.sum is not None:
new_tally._sum = np.zeros(new_tally.shape, dtype=np.float64)
new_tally._sum[diag_indices, :, :] = self.sum
if self.sum_sq is not None:
if not self.derived and self.sum_sq is not None:
new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64)
new_tally._sum_sq[diag_indices, :, :] = self.sum_sq
if self.mean is not None:

6
openmc/universe.py
View file

@ -1128,14 +1128,14 @@ class RectLattice(Lattice):
# For 2D Lattices
if len(self._dimension) == 2:
offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1]
offset = self._offsets[i[3]-1, i[2]-1, i[1]-1, distribcell_index-1]
offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path,
distribcell_index)
# For 3D Lattices
else:
offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1]
offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance(
offset = self._offsets[i[3]-1, i[2]-1, i[1]-1, distribcell_index-1]
offset += self._universes[i[3]-1][i[2]-1][i[1]-1].get_cell_instance(
path, distribcell_index)
return offset

2
setup.py
View file

@ -32,7 +32,7 @@ kwargs = {'name': 'openmc',
if have_setuptools:
kwargs.update({
# Required dependencies
'install_requires': ['numpy', 'h5py', 'matplotlib'],
'install_requires': ['numpy>=1.9', 'h5py', 'matplotlib'],
# Optional dependencies
'extras_require': {

770
src/ace.F90
View file

@ -1,19 +1,25 @@
module ace
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing, &
DistEnergy
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing
use constants
use endf, only: reaction_name, is_fission, is_disappearance
use error, only: fatal_error, warning
use fission, only: nu_total
use distribution_univariate, only: Uniform, Equiprobable, Tabular
use endf, only: is_fission, is_disappearance
use energy_distribution, only: TabularEquiprobable, LevelInelastic, &
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, NBodyPhaseSpace
use error, only: fatal_error, warning
use fission, only: nu_total
use global
use list_header, only: ListInt
use material_header, only: Material
use list_header, only: ListInt
use material_header, only: Material
use multipole, only: multipole_read
use multipole_header, only: max_L, max_poles, max_poly
use output, only: write_message
use set_header, only: SetChar
use string, only: to_str, to_lower
use output, only: write_message
use set_header, only: SetChar
use secondary_header, only: AngleEnergy
use secondary_correlated, only: CorrelatedAngleEnergy
use secondary_kalbach, only: KalbachMann
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str, to_lower
implicit none
@ -377,8 +383,8 @@ contains
else
call read_nu_data(nuc)
call read_reactions(nuc)
call read_angular_dist(nuc)
call read_energy_dist(nuc)
call read_angular_dist(nuc)
call read_unr_res(nuc)
end if
@ -584,9 +590,10 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
integer :: LAW
integer :: IDAT
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist
JXS2 = JXS(2)
JXS24 = JXS(24)
@ -729,11 +736,15 @@ contains
! Loop over all delayed neutron precursor groups
do i = 1, NPCR
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
LOCC = nint(XSS(LED + i - 1))
! Determine law and location of data
LAW = nint(XSS(LDIS + LOCC))
IDAT = nint(XSS(LDIS + LOCC + 1))
! read energy distribution data
edist => nuc % nu_d_edist(i)
call get_energy_dist(edist, LOCC, .true.)
call get_energy_dist(nuc%nu_d_edist(i)%obj, LAW, LDIS, IDAT, &
ZERO, ZERO)
end do
! =======================================================================
@ -806,8 +817,8 @@ contains
rxn%multiplicity = 1
rxn%threshold = 1
rxn%scatter_in_cm = .true.
rxn%has_angle_dist = .false.
rxn%has_energy_dist = .false.
allocate(rxn%secondary%distribution(1))
allocate(UncorrelatedAngleEnergy :: rxn%secondary%distribution(1)%obj)
end associate
! Add contribution of elastic scattering to total cross section
@ -822,10 +833,6 @@ contains
do i = 1, NMT
associate (rxn => nuc % reactions(i+1))
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
@ -955,86 +962,96 @@ contains
subroutine read_angular_dist(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
integer :: LOCB ! location of angular distribution for given MT
integer :: NE ! number of incoming energies
integer :: NP ! number of points for cosine distribution
integer :: LC ! locator
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
integer :: length ! length of data array to allocate
JXS8 = JXS(8)
JXS9 = JXS(9)
integer :: k ! index over energy distributions
integer :: interp
integer, allocatable :: LC(:) ! locator
! loop over all reactions with secondary neutrons -- NXS(5) does not include
! elastic scattering
do i = 1, NXS(5) + 1
associate (rxn => nuc%reactions(i))
! find location of angular distribution
LOCB = int(XSS(JXS8 + i - 1))
if (LOCB == -1) then
! Angular distribution data are specified through LAWi = 44 in the DLW
! block
cycle
elseif (LOCB == 0) then
! No angular distribution data are given for this reaction, isotropic
! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
cycle
end if
rxn % has_angle_dist = .true.
LOCB = int(XSS(JXS(8) + i - 1))
! allocate space for incoming energies and locations
NE = int(XSS(JXS9 + LOCB - 1))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
allocate(rxn % adist % location(NE))
! Angular distribution given as part of a correlated angle-energy distribution
if (LOCB == -1) cycle
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! No angular distribution data are given for this reaction, isotropic
! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
if (LOCB == 0) cycle
! determine dize of data block
length = 0
do j = 1, NE
LC = rxn % adist % location(j)
if (LC == 0) then
! isotropic
rxn % adist % type(j) = ANGLE_ISOTROPIC
elseif (LC > 0) then
! 32 equiprobable bins
rxn % adist % type(j) = ANGLE_32_EQUI
length = length + 33
elseif (LC < 0) then
! tabular distribution
rxn % adist % type(j) = ANGLE_TABULAR
NP = int(XSS(JXS9 + abs(LC)))
length = length + 2 + 3*NP
end if
end do
! Loop over each separate energy distribution. Even though there is only
! "one" angular distribution, it is repeated as many times as there are
! energy distributions for this reaction since the
! UncorrelatedAngleEnergy type holds one angle and energy distribution.
do k = 1, size(rxn%secondary%distribution)
select type (aedist => rxn%secondary%distribution(k)%obj)
type is (UncorrelatedAngleEnergy)
! allocate space for incoming energies and locations
NE = int(XSS(JXS(9) + LOCB - 1))
allocate(aedist%angle%energy(NE))
allocate(aedist%angle%distribution(NE))
allocate(LC(NE))
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! read incoming energy grid and location of nucs
XSS_index = JXS(9) + LOCB
aedist%angle%energy(:) = get_real(NE)
LC(:) = get_int(NE)
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
XSS_index = JXS9 + LOCB + 2 * NE
rxn % adist % data = get_real(length)
! determine dize of data block
do j = 1, NE
if (LC(j) == 0) then
! isotropic
allocate(Uniform :: aedist%angle%distribution(j)%obj)
select type (adist => aedist%angle%distribution(j)%obj)
type is (Uniform)
adist%a = -ONE
adist%b = ONE
end select
! change location pointers since they are currently relative to JXS(9)
LC = LOCB + 2 * NE + 1
do j = 1, NE
! For consistency, leave location as 0 if type is isotropic.
! This is not necessary for current correctness, but can avoid
! future issues
if (rxn % adist % location(j) /= 0) then
rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC
end if
elseif (LC(j) > 0) then
! 32 equiprobable bins
allocate(Equiprobable :: aedist%angle%distribution(j)%obj)
select type (adist => aedist%angle%distribution(j)%obj)
type is (Equiprobable)
allocate(adist%x(33))
end select
elseif (LC(j) < 0) then
! tabular distribution
allocate(Tabular :: aedist%angle%distribution(j)%obj)
end if
end do
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
do j = 1, NE
XSS_index = JXS(9) + abs(LC(j)) - 1
select type(adist => aedist%angle%distribution(j)%obj)
type is (Equiprobable)
adist%x(:) = get_real(33)
type is (Tabular)
! determine interpolation and number of points
interp = nint(XSS(XSS_index))
NP = nint(XSS(XSS_index + 1))
! Get probability density data
XSS_index = XSS_index + 2
allocate(adist%x(NP), adist%p(NP), adist%c(NP))
adist%x(:) = get_real(NP)
adist%p(:) = get_real(NP)
adist%c(:) = get_real(NP)
end select
end do
deallocate(LC)
end select
end do
end associate
end do
@ -1049,25 +1066,62 @@ contains
subroutine read_energy_dist(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: LED ! location of energy distribution locators
integer :: LOCC ! location of energy distributions for given MT
integer :: i ! loop index
LED = JXS(10)
integer :: n
integer :: IDAT ! locator for distribution data
integer :: LNW ! location of next energy law
integer :: LAW ! Type of energy law
! Loop over all reactions
do i = 1, NXS(5)
associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering
rxn % has_energy_dist = .true.
! Determine how many energy distributions are present for this reaction
LNW = nint(XSS(JXS(10) + i - 1))
n = 0
do while (LNW > 0)
n = n + 1
LNW = nint(XSS(JXS(11) + LNW - 1))
end do
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! Allocate space for distributions and probability of validity
associate (secondary => nuc%reactions(i + 1)%secondary)
allocate(secondary%applicability(n))
allocate(secondary%distribution(n))
! allocate energy distribution
allocate(rxn % edist)
LNW = nint(XSS(JXS(10) + i - 1))
n = 0
do while (LNW > 0)
n = n + 1
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
! Determine energy law and location of data
LAW = nint(XSS(JXS(11) + LNW))
IDAT = nint(XSS(JXS(11) + LNW + 1))
! Read probability of law validity
call secondary%applicability(n)%from_ace(XSS, JXS(11) + LNW + 2)
! Read energy law data
call get_energy_dist(secondary%distribution(n)%obj, LAW, &
JXS(11), IDAT, nuc%awr, nuc%reactions(i + 1)%Q_value)
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, when the angle/energy
! distribution was uncorrelated, no angle was actually sampled. With
! the refactor, an angle is always sampled for an uncorrelated
! distribution even when no angle distribution exists in the ACE file
! (isotropic is assumed). To preserve the RNG stream, we explicitly
! mark fission reactions so that we avoid the angle sampling.
if (any(nuc%reactions(i + 1)%MT == &
[N_FISSION, N_F, N_NF, N_2NF, N_3NF])) then
select type (aedist => secondary%distribution(n)%obj)
type is (UncorrelatedAngleEnergy)
aedist%fission = .true.
end select
end if
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
! Get locator for next distribution
LNW = nint(XSS(JXS(11) + LNW - 1))
end do
end associate
end do
@ -1079,281 +1133,319 @@ contains
! single reaction
!===============================================================================
recursive subroutine get_energy_dist(edist, loc_law, delayed_n)
type(DistEnergy), intent(inout) :: edist ! energy distribution
integer, intent(in) :: loc_law ! locator for data
logical, intent(in), optional :: delayed_n ! is this for delayed neutrons?
recursive subroutine get_energy_dist(aedist, law, LDIS, IDAT, awr, Q_value)
class(AngleEnergy), allocatable, intent(inout) :: aedist
integer, intent(in) :: law
integer, intent(in) :: LDIS
integer, intent(in) :: IDAT
real(8), intent(in) :: awr
real(8), intent(in) :: Q_value
integer :: LDIS ! location of all energy distributions
integer :: LNW ! location of next energy distribution if multiple
integer :: LAW ! secondary energy distribution law
integer :: i, j
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: IDAT ! location of first energy distribution for given MT
integer :: lc ! locator
integer :: length ! length of data to allocate
integer :: length_interp_data ! length of interpolation data
integer :: NP ! number of outgoing energies/angles
integer :: interp
integer, allocatable :: L(:) ! locations of distributions for each Ein
integer, allocatable :: LC(:) ! locations of distributions for each Ein
! determine location of energy distribution
if (present(delayed_n)) then
LDIS = JXS(27)
XSS_index = LDIS + IDAT - 1
if (law == 44) then
allocate(KalbachMann :: aedist)
elseif (law == 61) then
allocate(CorrelatedAngleEnergy :: aedist)
else
LDIS = JXS(11)
allocate(UncorrelatedAngleEnergy :: aedist)
end if
! locator for next law and information on this law
LNW = int(XSS(LDIS + loc_law - 1))
LAW = int(XSS(LDIS + loc_law))
IDAT = int(XSS(LDIS + loc_law + 1))
NR = int(XSS(LDIS + loc_law + 2))
edist % law = LAW
edist % p_valid % n_regions = NR
select type (aedist)
type is (UncorrelatedAngleEnergy)
! ========================================================================
! UNCORRELATED ENERGY DISTRIBUTIONS
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(edist % p_valid % nbt(NR))
allocate(edist % p_valid % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = LDIS + loc_law + 3
if (NR > 0) then
edist % p_valid % nbt = int(get_real(NR))
edist % p_valid % int = int(get_real(NR))
end if
! allocate space for law validity data
NE = int(XSS(LDIS + loc_law + 3 + 2*NR))
edist % p_valid % n_pairs = NE
allocate(edist % p_valid % x(NE))
allocate(edist % p_valid % y(NE))
length_interp_data = 5 + 2*(NR + NE)
! read law validity data
XSS_index = LDIS + loc_law + 4 + 2*NR
edist % p_valid % x = get_real(NE)
edist % p_valid % y = get_real(NE)
! Set index to beginning of IDAT array
lc = LDIS + IDAT - 2
! determine length of energy distribution
length = length_energy_dist(lc, LAW, loc_law, length_interp_data)
! allocate secondary energy distribution array
allocate(edist % data(length))
! read secondary energy distribution
XSS_index = lc + 1
edist % data = get_real(length)
! read next energy distribution if present
if (LNW > 0) then
allocate(edist % next)
call get_energy_dist(edist % next, LNW)
end if
end subroutine get_energy_dist
!===============================================================================
! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy
! distribution array based on the secondary energy law and location in XSS
!===============================================================================
function length_energy_dist(lc, law, LOCC, lid) result(length)
integer, intent(in) :: lc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
integer, intent(in) :: lid ! length of interpolation data
integer :: length ! length of energy distribution (LDAT)
integer :: i ! loop index for incoming energies
integer :: j ! loop index for outgoing energies
integer :: k ! dummy index in XSS
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
integer :: NP ! number of points in outgoing energy distribution
integer :: NMU ! number of points in outgoing cosine distribution
integer :: NRa ! number of interpolation regions for Watt 'a'
integer :: NEa ! number of energies for Watt 'a'
integer :: NRb ! number of interpolation regions for Watt 'b'
integer :: NEb ! number of energies for Watt 'b'
real(8), allocatable :: L(:) ! locations of distributions for each Ein
! initialize length
length = 0
select case (law)
case (1)
! Tabular equiprobable energy bins
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + NE))
length = 3 + 2*NR + NE + 3*NP*NE
case (2)
! Discrete photon energy
length = 2
case (3)
! Level scattering
length = 2
case (4)
! Continuous tabular distribution
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
select case (law)
case (1)
allocate(TabularEquiprobable :: aedist%energy)
select type (edist => aedist%energy)
type is (TabularEquiprobable)
NR = nint(XSS(XSS_index))
NE = nint(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for tabular equiprobable energy distributions.")
end if
end if
! determine length
NP = int(XSS(lc + length + 2))
length = length + 2 + 3*NP
edist%n_region = NR
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
! Read incoming energies for which outgoing energies are tabulated
allocate(edist%energy_in(NE))
XSS_index = XSS_index + 2 + 2*NR
edist%energy_in(:) = get_real(NE)
! Read outgoing energy tables
NP = nint(XSS(XSS_index))
allocate(edist%energy_out(NP, NE))
XSS_index = XSS_index + 1
do i = 1, NE
edist%energy_out(:, i) = get_real(NP)
end do
end select
case (3)
allocate(LevelInelastic :: aedist%energy)
select type (edist => aedist%energy)
type is (LevelInelastic)
edist%threshold = XSS(XSS_index)
edist%mass_ratio = XSS(XSS_index + 1)
end select
case (4)
allocate(ContinuousTabular :: aedist%energy)
select type (edist => aedist%energy)
type is (ContinuousTabular)
NR = nint(XSS(XSS_index))
XSS_index = XSS_index + 1
if (NR > 1) then
call fatal_error("Multiple interpolation regions not yet supported &
&for continuous tabular energy distributions.")
end if
edist%n_region = NR
! Read breakpoints and interpolation parameters
if (NR > 0) then
allocate(edist%breakpoints(NR))
allocate(edist%interpolation(NR))
edist%breakpoints(:) = get_int(NR)
edist%interpolation(:) = get_int(NR)
end if
! Read incoming energies for which outgoing energies are tabulated and
! locators
NE = nint(XSS(XSS_index))
XSS_index = XSS_index + 1
allocate(edist%energy_in(NE))
allocate(L(NE))
edist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Read outgoing energy tables
allocate(edist%energy_out(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
edist%energy_out(i)%interpolation = mod(interp, 10)
edist%energy_out(i)%n_discrete = (interp - &
edist%energy_out(i)%interpolation)/10
! check for discrete lines present
if (edist%energy_out(i)%n_discrete > 0) then
call fatal_error("Discrete lines in continuous tabular &
&distribution not yet supported")
end if
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(edist%energy_out(i)%e_out(NP))
allocate(edist%energy_out(i)%p(NP))
allocate(edist%energy_out(i)%c(NP))
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
edist%energy_out(i)%e_out(:) = get_real(NP)
edist%energy_out(i)%p(:) = get_real(NP)
edist%energy_out(i)%c(:) = get_real(NP)
end do
deallocate(L)
end select
case (7)
allocate(MaxwellEnergy :: aedist%energy)
select type (edist => aedist%energy)
type is (MaxwellEnergy)
call edist%theta%from_ace(XSS, XSS_index)
edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
2*edist%theta%n_pairs)
end select
case (9)
allocate(Evaporation :: aedist%energy)
select type(edist => aedist%energy)
type is (Evaporation)
call edist%theta%from_ace(XSS, XSS_index)
edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
2*edist%theta%n_pairs)
end select
case (11)
allocate(WattEnergy :: aedist%energy)
select type(edist => aedist%energy)
type is (WattEnergy)
call edist%a%from_ace(XSS, XSS_index)
XSS_index = XSS_index + 2 + 2*edist%a%n_regions + 2*edist%a%n_pairs
call edist%b%from_ace(XSS, XSS_index)
XSS_index = XSS_index + 2 + 2*edist%b%n_regions + 2*edist%b%n_pairs
edist%u = XSS(XSS_index)
end select
case (66)
allocate(NBodyPhaseSpace :: aedist%energy)
select type(edist => aedist%energy)
type is (NBodyPhaseSpace)
edist%n_bodies = int(XSS(XSS_index))
edist%mass_ratio = XSS(XSS_index + 1)
edist%A = awr
edist%Q = Q_value
end select
end select
type is (KalbachMann)
! ========================================================================
! CORRELATED KALBACH-MANN DISTRIBUTION
NR = int(XSS(XSS_index))
NE = int(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for Kalbach-Mann energy distributions.")
end if
aedist%n_region = NR
! Read incoming energies for which outgoing energies are tabulated and locators
allocate(aedist%energy_in(NE))
allocate(L(NE))
XSS_index = XSS_index + 2 + 2*NR
aedist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Read outgoing energy tables
allocate(aedist%table(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
aedist%table(i)%interpolation = mod(interp, 10)
aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
! check for discrete lines present
if (aedist%table(i)%n_discrete > 0) then
call fatal_error("Discrete lines in Kalbach-Mann distribution not &
&yet supported")
end if
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(aedist%table(i)%e_out(NP))
allocate(aedist%table(i)%p(NP))
allocate(aedist%table(i)%c(NP))
allocate(aedist%table(i)%r(NP))
allocate(aedist%table(i)%a(NP))
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
aedist%table(i)%e_out(:) = get_real(NP)
aedist%table(i)%p(:) = get_real(NP)
aedist%table(i)%c(:) = get_real(NP)
aedist%table(i)%r(:) = get_real(NP)
aedist%table(i)%a(:) = get_real(NP)
end do
deallocate(L)
case (5)
! General evaporation spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
NP = int(XSS(lc + 3 + 2*NR + 2*NE))
length = 3 + 2*NR + 2*NE + NP
type is (CorrelatedAngleEnergy)
! ========================================================================
! CORRELATED ANGLE-ENERGY DISTRIBUTION
case (7)
! Maxwell fission spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
NR = int(XSS(XSS_index))
NE = int(XSS(XSS_index + 1 + 2*NR))
if (NR > 0) then
call fatal_error("Multiple interpolation regions not yet supported &
&for correlated angle-energy distributions.")
end if
aedist%n_region = NR
case (9)
! Evaporation spectrum
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
length = 3 + 2*NR + 2*NE
case (11)
! Watt spectrum
NRa = int(XSS(lc + 1))
NEa = int(XSS(lc + 2 + 2*NRa))
NRb = int(XSS(lc + 3 + 2*(NRa+NEa)))
NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb)))
length = 5 + 2*(NRa + NEa + NRb + NEb)
case (44)
! Kalbach-Mann correlated scattering
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
! Read incoming energies for which outgoing energies are tabulated and
! locators
allocate(aedist%energy_in(NE))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
XSS_index = XSS_index + 2 + 2*NR
aedist%energy_in(:) = get_real(NE)
L(:) = get_int(NE)
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
end if
end if
NP = int(XSS(lc + length + 2))
length = length + 2 + 5*NP
! Read outgoing energy tables
allocate(aedist%table(NE))
do i = 1, NE
! Determine interpolation and number of discrete points
XSS_index = LDIS + L(i) - 1
interp = nint(XSS(XSS_index))
aedist%table(i)%interpolation = mod(interp, 10)
aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
! adjust location for this block
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
deallocate(L)
case (61)
! Correlated energy and angle distribution
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
do i = 1,NE
! Some older data sets use the same LDAT for multiple Ein tables.
! If this is the case, we should skip incrementing length when it is
! not needed.
if (i < NE) then
if (any(L(i) == L(i + 1: NE))) then
! adjust locators for energy distribution
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
cycle
end if
! check for discrete lines present
if (aedist%table(i)%n_discrete > 0) then
call fatal_error("Discrete lines in correlated angle-energy &
&distribution not yet supported")
end if
! outgoing energy distribution
NP = int(XSS(lc + length + 2))
! Determine number of points and allocate space
NP = nint(XSS(XSS_index + 1))
allocate(aedist%table(i)%e_out(NP))
allocate(aedist%table(i)%p(NP))
allocate(aedist%table(i)%c(NP))
allocate(LC(NP))
! adjust locators for angular distribution
! Read tabular PDF for outgoing energy
XSS_index = XSS_index + 2
aedist%table(i)%e_out(:) = get_real(NP)
aedist%table(i)%p(:) = get_real(NP)
aedist%table(i)%c(:) = get_real(NP)
LC(:) = get_int(NP)
! allocate angular distributions for each incoming/outgoing energy
allocate(aedist%table(i)%angle(NP))
do j = 1, NP
k = lc + length + 2 + 3*NP + j
if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
if (LC(j) == 0) then
! isotropic
allocate(Uniform :: aedist%table(i)%angle(j)%obj)
select type (adist => aedist%table(i)%angle(j)%obj)
type is (Uniform)
adist%a = -ONE
adist%b = ONE
end select
elseif (LC(j) > 0) then
! tabular distribution
allocate(Tabular :: aedist%table(i)%angle(j)%obj)
end if
end do
length = length + 2 + 4*NP
! read angular distributions
do j = 1, NP
! outgoing angle distribution -- NMU here is actually
! referred to as NP in the MCNP documentation
NMU = int(XSS(lc + length + 2))
length = length + 2 + 3*NMU
XSS_index = LDIS + abs(LC(j)) - 1
select type(adist => aedist%table(i)%angle(j)%obj)
type is (Tabular)
! determine interpolation and number of points
interp = nint(XSS(XSS_index))
NP = nint(XSS(XSS_index + 1))
! Get probability density data
XSS_index = XSS_index + 2
allocate(adist%x(NP), adist%p(NP), adist%c(NP))
adist%x(:) = get_real(NP)
adist%p(:) = get_real(NP)
adist%c(:) = get_real(NP)
end select
end do
deallocate(LC)
! adjust locators for energy distribution
j = lc + 2 + 2*NR + NE + i
XSS(j) = XSS(j) - LOCC - lid
end do
deallocate(L)
case (66)
! N-body phase space distribution
length = 2
case (67)
! Laboratory energy-angle law
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
! Before progressing, check to see if data set uses L(I) values
! in a way inconsistent with the current form of the ACE Format Guide
! (MCNP5 Manual, Vol 3)
allocate(L(NE))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Don't currently do anything with L
deallocate(L)
! Continue with finding data length
NMU = int(XSS(lc + 4 + 2*NR + 2*NE))
length = 4 + 2*(NR + NE + NMU)
end select
end function length_energy_dist
end subroutine get_energy_dist
!===============================================================================
! READ_UNR_RES reads in unresolved resonance probability tables if present.

85
src/ace_header.F90
View file

@ -1,45 +1,14 @@
module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use multipole_header, only: MultipoleArray
use stl_vector, only: VectorInt
use secondary_header, only: SecondaryDistribution, AngleEnergyContainer
use stl_vector, only: VectorInt
implicit none
!===============================================================================
! DISTANGLE contains data for a tabular secondary angle distribution whether it
! be tabular or 32 equiprobable cosine bins
!===============================================================================
type DistAngle
integer :: n_energy ! # of incoming energies
real(8), allocatable :: energy(:) ! incoming energy grid
integer, allocatable :: type(:) ! type of distribution
integer, allocatable :: location(:) ! location of each table
real(8), allocatable :: data(:) ! angular distribution data
end type DistAngle
!===============================================================================
! DISTENERGY contains data for a secondary energy distribution for all
! scattering laws
!===============================================================================
type DistEnergy
integer :: law ! secondary distribution law
type(Tab1) :: p_valid ! probability of law validity
real(8), allocatable :: data(:) ! energy distribution data
! For reactions that may have multiple energy distributions such as (n,2n),
! this pointer allows multiple laws to be stored
type(DistEnergy), pointer :: next => null()
! Type-Bound procedures
contains
procedure :: clear => distenergy_clear ! Deallocates DistEnergy
end type DistEnergy
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
@ -54,10 +23,7 @@ module ace_header
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
real(8), allocatable :: sigma(:) ! Cross section values
logical :: has_angle_dist ! Angle distribution present?
logical :: has_energy_dist ! Energy distribution present?
type(DistAngle) :: adist ! Secondary angular distribution
type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution
type(SecondaryDistribution) :: secondary
! Type-Bound procedures
contains
@ -138,7 +104,7 @@ module ace_header
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(DistEnergy), pointer :: nu_d_edist(:) => null()
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
! Unresolved resonance data
logical :: urr_present
@ -292,37 +258,14 @@ module ace_header
contains
!===============================================================================
! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
!===============================================================================
recursive subroutine distenergy_clear(this)
class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
if (associated(this % next)) then
! recursively clear this item
call this % next % clear()
deallocate(this % next)
end if
end subroutine distenergy_clear
!===============================================================================
! REACTION_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine reaction_clear(this)
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)
end if
end subroutine reaction_clear
!===============================================================================
@ -330,21 +273,11 @@ module ace_header
!===============================================================================
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
class(Nuclide), intent(inout) :: this
integer :: i ! Loop counter
if (associated(this % nu_d_edist)) then
do i = 1, size(this % nu_d_edist)
call this % nu_d_edist(i) % clear()
end do
deallocate(this % nu_d_edist)
end if
if (associated(this % urr_data)) then
deallocate(this % urr_data)
end if
if (associated(this % urr_data)) deallocate(this % urr_data)
if (this % mp_present) then
deallocate(this % multipole)

63
src/angle_distribution.F90 Normal file
View file

@ -0,0 +1,63 @@
module angle_distribution
use constants, only: ZERO, ONE
use distribution_univariate, only: DistributionContainer
use random_lcg, only: prn
use search, only: binary_search
implicit none
private
!===============================================================================
! ANGLEDISTRIBUTION represents an angular distribution that is to be used in an
! uncorrelated angle-energy distribution. This occurs whenever the angle
! distrbution is given in File 4 in an ENDF file. The distribution of angles
! depends on the incoming energy of the neutron, so this type stores a
! distribution for each of a set of incoming energies.
!===============================================================================
type, public :: AngleDistribution
real(8), allocatable :: energy(:)
type(DistributionContainer), allocatable :: distribution(:)
contains
procedure :: sample => angle_sample
end type AngleDistribution
contains
function angle_sample(this, E) result(mu)
class(AngleDistribution), intent(in) :: this
real(8), intent(in) :: E ! incoming energy
real(8) :: mu ! sampled cosine of scattering angle
integer :: i ! index on incoming energy grid
integer :: n ! number of incoming energies
real(8) :: r ! interpolation factor on incoming energy grid
! Determine number of incoming energies
n = size(this%energy)
! Find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
if (E < this%energy(1)) then
i = 1
r = ZERO
elseif (E > this%energy(n)) then
i = n - 1
r = ONE
else
i = binary_search(this%energy, n, E)
r = (E - this%energy(i))/(this%energy(i+1) - this%energy(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) i = i + 1
! Sample i-th distribution
mu = this%distribution(i)%obj%sample()
! Make sure mu is in range [-1,1]
if (abs(mu) > ONE) mu = sign(ONE, mu)
end function angle_sample
end module angle_distribution

29
src/distribution_univariate.F90
View file

@ -1,7 +1,7 @@
module distribution_univariate
use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR, MAX_LINE_LEN, &
MAX_WORD_LEN
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
@ -78,6 +78,12 @@ module distribution_univariate
procedure :: initialize => tabular_initialize
end type Tabular
type, extends(Distribution) :: Equiprobable
real(8), allocatable :: x(:)
contains
procedure :: sample => equiprobable_sample
end type Equiprobable
contains
function discrete_sample(this) result(x)
@ -238,6 +244,25 @@ contains
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
function equiprobable_sample(this) result(x)
class(Equiprobable), intent(in) :: this
real(8) :: x
integer :: i
integer :: n
real(8) :: r
real(8) :: xl, xr
n = size(this%x)
r = prn()
i = 1 + int((n - 1)*r)
xl = this%x(i)
xr = this%x(i+1)
x = xl + ((n - 1)*r - i + ONE) * (xr - xl)
end function equiprobable_sample
subroutine distribution_from_xml(dist, node_dist)
class(Distribution), allocatable, intent(inout) :: dist
type(Node), pointer :: node_dist

47
src/endf.F90
View file

@ -17,6 +17,53 @@ contains
character(20) :: string
select case (MT)
! Special reactions for tallies
case (SCORE_FLUX)
string = "flux"
case (SCORE_TOTAL)
string = "total"
case (SCORE_SCATTER)
string = "scatter"
case (SCORE_NU_SCATTER)
string = "nu-scatter"
case (SCORE_SCATTER_N)
string = "scatter-n"
case (SCORE_SCATTER_PN)
string = "scatter-pn"
case (SCORE_NU_SCATTER_N)
string = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
string = "nu-scatter-pn"
case (SCORE_TRANSPORT)
string = "transport"
case (SCORE_N_1N)
string = "n1n"
case (SCORE_ABSORPTION)
string = "absorption"
case (SCORE_FISSION)
string = "fission"
case (SCORE_NU_FISSION)
string = "nu-fission"
case (SCORE_DELAYED_NU_FISSION)
string = "delayed-nu-fission"
case (SCORE_KAPPA_FISSION)
string = "kappa-fission"
case (SCORE_CURRENT)
string = "current"
case (SCORE_FLUX_YN)
string = "flux-yn"
case (SCORE_TOTAL_YN)
string = "total-yn"
case (SCORE_SCATTER_YN)
string = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
string = "nu-scatter-yn"
case (SCORE_EVENTS)
string = "events"
case (SCORE_INVERSE_VELOCITY)
string = "inverse-velocity"
! Normal ENDF-based reactions
case (TOTAL_XS)
string = '(n,total)'
case (ELASTIC)

34
src/endf_header.F90
View file

@ -13,6 +13,40 @@ module endf_header
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
contains
procedure :: from_ace
end type Tab1
contains
subroutine from_ace(this, xss, idx)
class(Tab1), intent(inout) :: this
real(8), intent(in) :: xss(:)
integer, intent(in) :: idx
integer :: nr, ne
! Determine number of regions
nr = nint(xss(idx))
this%n_regions = nr
! Read interpolation region data
if (nr > 0) then
allocate(this%nbt(nr))
allocate(this%int(nr))
this%nbt(:) = nint(xss(idx + 1 : idx + nr))
this%int(:) = nint(xss(idx + nr + 1 : idx + 2*nr))
end if
! Determine number of pairs
ne = int(XSS(idx + 2*nr + 1))
this%n_pairs = ne
! Read (x,y) pairs
allocate(this%x(ne))
allocate(this%y(ne))
this%x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne)
this%y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne)
end subroutine from_ace
end module endf_header

429
src/energy_distribution.F90 Normal file
View file

@ -0,0 +1,429 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a
! function of the incident energy of a projectile. Each derived type must
! implement a sample() function that returns a sampled outgoing energy given an
! incoming energy
!===============================================================================
type, abstract :: EnergyDistribution
contains
procedure(iSampleEnergy), deferred :: sample
end type EnergyDistribution
abstract interface
function iSampleEnergy(this, E_in) result(E_out)
import EnergyDistribution
class(EnergyDistribution), intent(in) :: this
real(8), intent(in) :: E_in
real(8) :: E_out
end function iSampleEnergy
end interface
type :: EnergyDistributionContainer
class(EnergyDistribution), allocatable :: obj
end type EnergyDistributionContainer
!===============================================================================
! Derived classes
!===============================================================================
!===============================================================================
! TABULAREQUIPROBABLE represents an energy distribution with tabular
! equiprobable energy bins as given in ACE law 1. This is an older
! representation that has largely been replaced with ACE laws 4, 44, and 61.
!===============================================================================
type, extends(EnergyDistribution) :: TabularEquiprobable
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
real(8), allocatable :: energy_out(:,:) ! table of outgoing energies for
! each incoming energy
contains
procedure :: sample => equiprobable_sample
end type TabularEquiprobable
!===============================================================================
! LEVELINELASTIC gives the energy distribution for level inelastic scattering by
! neutrons as in ENDF MT=51--90.
!===============================================================================
type, extends(EnergyDistribution) :: LevelInelastic
real(8) :: threshold
real(8) :: mass_ratio
contains
procedure :: sample => level_inelastic_sample
end type LevelInelastic
!===============================================================================
! CONTINUOUSTABULAR gives an energy distribution represented as a tabular
! distribution with histogram or linear-linear interpolation. This corresponds
! to ACE law 4, which NJOY produces for a number of ENDF energy distributions.
!===============================================================================
type CTTable
integer :: interpolation
integer :: n_discrete
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
end type CTTable
type, extends(EnergyDistribution) :: ContinuousTabular
integer :: n_region
integer, allocatable :: breakpoints(:)
integer, allocatable :: interpolation(:)
real(8), allocatable :: energy_in(:)
type(CTTable), allocatable :: energy_out(:)
contains
procedure :: sample => continuous_sample
end type ContinuousTabular
!===============================================================================
! MAXWELLENERGY gives the energy distribution of neutrons emitted from a Maxwell
! fission spectrum. This corresponds to ACE law 7 and ENDF File 5, LF=7.
!===============================================================================
type, extends(EnergyDistribution) :: MaxwellEnergy
type(Tab1) :: theta ! incoming-energy-dependent parameter
real(8) :: u ! restriction energy
contains
procedure :: sample => maxwellenergy_sample
end type MaxwellEnergy
!===============================================================================
! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and
! ENDF File 5, LF=9.
!===============================================================================
type, extends(EnergyDistribution) :: Evaporation
type(Tab1) :: theta
real(8) :: u
contains
procedure :: sample => evaporation_sample
end type Evaporation
!===============================================================================
! WATTENERGY gives the energy distribution of neutrons emitted from a Watt
! fission spectrum. This corresponds to ACE law 11 and ENDF File 5, LF=11.
!===============================================================================
type, extends(EnergyDistribution) :: WattEnergy
type(Tab1) :: a
type(Tab1) :: b
real(8) :: u
contains
procedure :: sample => watt_sample
end type WattEnergy
!===============================================================================
! NBODYPHASESPACE gives the energy distribution for particles emitted from
! neutron and charged-particle reactions. This corresponds to ACE law 66 and
! ENDF File 6, LAW=6.
!===============================================================================
type, extends(EnergyDistribution) :: NBodyPhaseSpace
integer :: n_bodies
real(8) :: mass_ratio
real(8) :: A
real(8) :: Q
contains
procedure :: sample => nbody_sample
end type NBodyPhaseSpace
contains
function equiprobable_sample(this, E_in) result(E_out)
class(TabularEquiprobable), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
! Determine number of incoming/outgoing energies
n_energy_in = size(this%energy_in)
n_energy_out = size(this%energy_out, 1)
! Determine index on incoming energy grid and interpolation factor
i = binary_search(this%energy_in, size(this%energy_in), E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
! Sample outgoing energy bin
k = 1 + int(n_energy_out * prn())
! Determine E_1 and E_K
E_i_1 = this%energy_out(1, i)
E_i_K = this%energy_out(n_energy_out, i)
E_i1_1 = this%energy_out(1, i+1)
E_i1_K = this%energy_out(n_energy_out, i+1)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! Randomly select between the outgoing table for incoming energy E_i and
! E_(i+1)
if (prn() < r) then
l = i + 1
else
l = i
end if
! Determine E_l_k and E_l_k+1
E_l_k = this%energy_out(k, l)
E_l_k1 = this%energy_out(k+1, l)
! Determine E' (denoted here as E_out)
E_out = E_l_k + prn()*(E_l_k1 - E_l_k)
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
end function equiprobable_sample
function level_inelastic_sample(this, E_in) result(E_out)
class(LevelInelastic), intent(in) :: this
real(8), intent(in) :: E_in
real(8) :: E_out
E_out = this%mass_ratio*(E_in - this%threshold)
end function level_inelastic_sample
function continuous_sample(this, E_in) result(E_out)
class(ContinuousTabular), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
logical :: histogram_interp ! whether histogram interpolation is used
! Read number of interpolation regions and incoming energies
if (this%n_region == 1) then
histogram_interp = (this%interpolation(1) == 1)
else
histogram_interp = .false.
end if
! Find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (histogram_interp) then
l = i
else
if (r > prn()) then
l = i + 1
else
l = i
end if
end if
! Interpolation for energy E1 and EK
n_energy_out = size(this%energy_out(i)%e_out)
E_i_1 = this%energy_out(i)%e_out(1)
E_i_K = this%energy_out(i)%e_out(n_energy_out)
n_energy_out = size(this%energy_out(i+1)%e_out)
E_i1_1 = this%energy_out(i+1)%e_out(1)
E_i1_K = this%energy_out(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! Determine outgoing energy bin
n_energy_out = size(this%energy_out(l)%e_out)
r1 = prn()
c_k = this%energy_out(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%energy_out(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! Check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%energy_out(l)%e_out(k)
p_l_k = this%energy_out(l)%p(k)
if (this%energy_out(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%energy_out(l)%e_out(k+1)
p_l_k1 = this%energy_out(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
end if
! Now interpolate between incident energy bins i and i + 1
if (.not. histogram_interp) then
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
end if
end function continuous_sample
function maxwellenergy_sample(this, E_in) result(E_out)
class(MaxwellEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! Maxwell distribution parameter
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
do
! Sample maxwell fission spectrum
E_out = maxwell_spectrum(theta)
! Accept energy based on restriction energy
if (E_out <= E_in - this%u) exit
end do
end function maxwellenergy_sample
function evaporation_sample(this, E_in) result(E_out)
class(Evaporation), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: theta ! evaporation spectrum parameter
real(8) :: x, y, v
! Get temperature corresponding to incoming energy
theta = interpolate_tab1(this%theta, E_in)
y = (E_in - this%U)/theta
v = 1 - exp(-y)
! Sample outgoing energy based on evaporation spectrum probability
! density function
do
x = -log((ONE - v*prn())*(ONE - v*prn()))
if (x <= y) exit
end do
E_out = x*theta
end function evaporation_sample
function watt_sample(this, E_in) result(E_out)
class(WattEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: a, b ! Watt spectrum parameters
! Determine Watt parameter 'a' from tabulated function
a = interpolate_tab1(this%a, E_in)
! Determine Watt parameter 'b' from tabulated function
b = interpolate_tab1(this%b, E_in)
do
! Sample energy-dependent Watt fission spectrum
E_out = watt_spectrum(a, b)
! Accept energy based on restriction energy
if (E_out <= E_in - this%u) exit
end do
end function watt_sample
function nbody_sample(this, E_in) result(E_out)
class(NBodyPhaseSpace), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8) :: E_out ! sampled outgoing energy
real(8) :: Ap ! total mass of particles in neutron masses
real(8) :: E_max ! maximum possible COM energy
real(8) :: x, y, v
real(8) :: r1, r2, r3, r4, r5, r6
! Determine E_max parameter
Ap = this%mass_ratio
E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q)
! x is essentially a Maxwellian distribution
x = maxwell_spectrum(ONE)
select case (this%n_bodies)
case (3)
y = maxwell_spectrum(ONE)
case (4)
r1 = prn()
r2 = prn()
r3 = prn()
y = -log(r1*r2*r3)
case (5)
r1 = prn()
r2 = prn()
r3 = prn()
r4 = prn()
r5 = prn()
r6 = prn()
y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
end select
! Now determine v and E_out
v = x/(x+y)
E_out = E_max * v
end function nbody_sample
end module energy_distribution

1
src/global.F90
View file

@ -22,7 +22,6 @@ module global
#endif
implicit none
save
! ============================================================================
! GEOMETRY-RELATED VARIABLES

29
src/input_xml.F90
View file

@ -5,6 +5,7 @@ module input_xml
use dict_header, only: DictIntInt, ElemKeyValueCI
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
@ -3281,7 +3282,7 @@ contains
call fatal_error("Cannot tally absorption rate with an outgoing &
&energy filter.")
end if
case ('fission')
case ('fission', '18')
t % score_bins(j) = SCORE_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally fission rate with an outgoing &
@ -3461,6 +3462,32 @@ contains
! Deallocate temporary string array of scores
deallocate(sarray)
! Check that no duplicate scores exist
j = 1
do while (j < n_scores)
! Determine number of bins for scores with expansions
n_order = t % moment_order(j)
select case (t % score_bins(j))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
n_bins = n_order + 1
case (SCORE_FLUX_YN, SCORE_TOTAL_YN, SCORE_SCATTER_YN, &
SCORE_NU_SCATTER_YN)
n_bins = (n_order + 1)**2
case default
n_bins = 1
end select
do k = j + n_bins, n_scores
if (t % score_bins(j) == t % score_bins(k) .and. &
t % moment_order(j) == t % moment_order(k)) then
call fatal_error("Duplicate score of type '" // trim(&
reaction_name(t % score_bins(j))) // "' found in tally " &
// trim(to_str(t % id)))
end if
end do
j = j + n_bins
end do
else
call fatal_error("No <scores> specified on tally " &
&// trim(to_str(t % id)) // ".")

1
src/interpolation.F90
View file

@ -2,7 +2,6 @@ module interpolation
use constants
use endf_header, only: Tab1
use error, only: fatal_error
use search, only: binary_search
use string, only: to_str

40
src/output.F90
View file

@ -322,14 +322,8 @@ contains
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_total ! memory used by nuclide (bytes)
integer :: size_angle_total ! total memory used for angle dist. (bytes)
integer :: size_energy_total ! total memory used for energy dist. (bytes)
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_angle ! memory used for an angle distribution (bytes)
integer :: size_energy ! memory used for a energy distributions (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
character(11) :: law ! secondary energy distribution law
type(UrrData), pointer :: urr
! set default unit for writing information
@ -340,8 +334,6 @@ contains
end if
! Initialize totals
size_angle_total = 0
size_energy_total = 0
size_urr = 0
size_xs = 0
@ -356,33 +348,15 @@ contains
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, nuc % n_reaction
associate (rxn => nuc % reactions(i))
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution and law
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
law = to_str(rxn % edist % law)
else
size_energy = 0
law = 'None'
end if
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
rxn % threshold
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
size_angle_total = size_angle_total + size_angle
size_energy_total = size_energy_total + size_energy
end associate
end do
@ -408,19 +382,11 @@ contains
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end if
! Calculate total memory
size_total = size_xs + size_angle_total + size_energy_total + size_urr
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Secondary angle distributions = ' // &
trim(to_str(size_angle_total)) // ' bytes'
write(unit_,*) ' Secondary energy distributions = ' // &
trim(to_str(size_energy_total)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)

2
src/particle_header.F90
View file

@ -89,7 +89,7 @@ module particle_header
logical :: write_track = .false.
! Secondary particles created
integer :: n_secondary = 0
integer(8) :: n_secondary = 0
type(Bank) :: secondary_bank(MAX_SECONDARY)
contains

1015
src/physics.F90
View file

File diff suppressed because it is too large Load diff

4
src/relaxng/tallies.rnc
View file

@ -19,8 +19,8 @@ element tallies {
(element id { xsd:int } | attribute id { xsd:int }) &
(element name { xsd:string { maxLength="52" } } |
attribute name { xsd:string { maxLength="52" } })? &
(element estimator { ( "analog" | "tracklength" ) } |
attribute estimator { ( "analog" | "tracklength" ) })? &
(element estimator { ( "analog" | "tracklength" | "collision" ) } |
attribute estimator { ( "analog" | "tracklength" | "collision" ) })? &
element filter {
(element type { ( "cell" | "cellborn" | "material" | "universe" |
"surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" |

2
src/relaxng/tallies.rng
View file

@ -120,12 +120,14 @@
<choice>
<value>analog</value>
<value>tracklength</value>
<value>collision</value>
</choice>
</element>
<attribute name="estimator">
<choice>
<value>analog</value>
<value>tracklength</value>
<value>collision</value>
</choice>
</attribute>
</choice>

1
src/search.F90
View file

@ -1,7 +1,6 @@
module search
use constants
use error, only: fatal_error
implicit none

148
src/secondary_correlated.F90 Normal file
View file

@ -0,0 +1,148 @@
module secondary_correlated
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This
! corresponds to ACE law 61 and ENDF File 6, LAW=1, LANG/=2.
!===============================================================================
type AngleEnergyTable
integer :: interpolation
integer :: n_discrete
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
type(DistributionContainer), allocatable :: angle(:)
end type AngleEnergyTable
type, extends(AngleEnergy) :: CorrelatedAngleEnergy
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(AngleEnergyTable), allocatable :: table(:) ! outgoing E/mu distributions
contains
procedure :: sample => correlated_sample
end type CorrelatedAngleEnergy
contains
subroutine correlated_sample(this, E_in, E_out, mu)
class(CorrelatedAngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sapmled scattering cosine
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, an isotropic polar cosine was
! always sampled but then overwritten with the polar cosine sampled from the
! correlated distribution. To preserve the random number stream, we keep
! this dummy sampling here but can remove it later (will change answers)
mu = TWO*prn() - ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) then
l = i + 1
else
l = i
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
! Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1) then
mu = this%table(l)%angle(k)%obj%sample()
else
mu = this%table(l)%angle(k + 1)%obj%sample()
end if
end subroutine correlated_sample
end module secondary_correlated

78
src/secondary_header.F90 Normal file
View file

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

164
src/secondary_kalbach.F90 Normal file
View file

@ -0,0 +1,164 @@
module secondary_kalbach
use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! KalbachMann represents a correlated angle-energy distribution with the angular
! distribution represented using Kalbach-Mann systematics. This corresponds to
! ACE law 44 and ENDF File 6, LAW=1, LANG=2.
!===============================================================================
type KalbachMannTable
integer :: n_discrete
integer :: interpolation
real(8), allocatable :: e_out(:)
real(8), allocatable :: p(:)
real(8), allocatable :: c(:)
real(8), allocatable :: r(:)
real(8), allocatable :: a(:)
end type KalbachMannTable
type, extends(AngleEnergy) :: KalbachMann
integer :: n_region ! number of interpolation regions
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
integer, allocatable :: interpolation(:) ! interpolation region codes
real(8), allocatable :: energy_in(:) ! incoming energies
type(KalbachMannTable), allocatable :: table(:) ! outgoing E/mu parameters
contains
procedure :: sample => kalbachmann_sample
end type KalbachMann
contains
subroutine kalbachmann_sample(this, E_in, E_out, mu)
class(KalbachMann), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
integer :: i, k, l ! indices
integer :: n_energy_in ! number of incoming energies
integer :: n_energy_out ! number of outgoing energies
real(8) :: r ! interpolation factor on incoming energy
real(8) :: r1 ! random number on [0,1)
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
real(8) :: c_k, c_k1 ! cumulative probability
real(8) :: km_r, km_a ! Kalbach-Mann parameters
real(8) :: T
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, an isotropic polar cosine was
! always sampled but then overwritten with the polar cosine sampled from the
! correlated distribution. To preserve the random number stream, we keep
! this dummy sampling here but can remove it later (will change answers)
mu = TWO*prn() - ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! find energy bin and calculate interpolation factor -- if the energy is
! outside the range of the tabulated energies, choose the first or last bins
n_energy_in = size(this%energy_in)
if (E_in < this%energy_in(1)) then
i = 1
r = ZERO
elseif (E_in > this%energy_in(n_energy_in)) then
i = n_energy_in - 1
r = ONE
else
i = binary_search(this%energy_in, n_energy_in, E_in)
r = (E_in - this%energy_in(i)) / &
(this%energy_in(i+1) - this%energy_in(i))
end if
! Sample between the ith and (i+1)th bin
if (r > prn()) then
l = i + 1
else
l = i
end if
! interpolation for energy E1 and EK
n_energy_out = size(this%table(i)%e_out)
E_i_1 = this%table(i)%e_out(1)
E_i_K = this%table(i)%e_out(n_energy_out)
n_energy_out = size(this%table(i+1)%e_out)
E_i1_1 = this%table(i+1)%e_out(1)
E_i1_K = this%table(i+1)%e_out(n_energy_out)
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
E_K = E_i_K + r*(E_i1_K - E_i_K)
! determine outgoing energy bin
n_energy_out = size(this%table(l)%e_out)
r1 = prn()
c_k = this%table(l)%c(1)
do k = 1, n_energy_out - 1
c_k1 = this%table(l)%c(k+1)
if (r1 < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, n_energy_out - 1)
E_l_k = this%table(l)%e_out(k)
p_l_k = this%table(l)%p(k)
if (this%table(l)%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_l_k > ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k)
km_a = this%table(l)%a(k)
elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
! Linear-linear interpolation
E_l_k1 = this%table(l)%e_out(k+1)
p_l_k1 = this%table(l)%p(k+1)
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
if (frac == ZERO) then
E_out = E_l_k + (r1 - c_k)/p_l_k
else
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
TWO*frac*(r1 - c_k))) - p_l_k)/frac
end if
! Determine Kalbach-Mann parameters
km_r = this%table(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%r(k+1) - this%table(l)%r(k))
km_a = this%table(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
(this%table(l)%a(k+1) - this%table(l)%a(k))
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
else
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
end if
! Sampled correlated angle from Kalbach-Mann parameters
if (prn() > km_r) then
T = (TWO*prn() - ONE) * sinh(km_a)
mu = log(T + sqrt(T*T + ONE))/km_a
else
r1 = prn()
mu = log(r1*exp(km_a) + (ONE - r1)*exp(-km_a))/km_a
end if
end subroutine kalbachmann_sample
end module secondary_kalbach

48
src/secondary_uncorrelated.F90 Normal file
View file

@ -0,0 +1,48 @@
module secondary_uncorrelated
use angle_distribution, only: AngleDistribution
use constants, only: ONE, TWO
use energy_distribution, only: EnergyDistribution
use secondary_header, only: AngleEnergy
use random_lcg, only: prn
!===============================================================================
! UNCORRELATEDANGLEENERGY represents an uncorrelated angle-energy
! distribution. This corresponds to when an energy distribution is given in ENDF
! File 5/6 and an angular distribution is given in ENDF File 4.
!===============================================================================
type, extends(AngleEnergy) :: UncorrelatedAngleEnergy
logical :: fission = .false.
type(AngleDistribution) :: angle
class(EnergyDistribution), allocatable :: energy
contains
procedure :: sample => uncorrelated_sample
end type UncorrelatedAngleEnergy
contains
subroutine uncorrelated_sample(this, E_in, E_out, mu)
class(UncorrelatedAngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in ! incoming energy
real(8), intent(out) :: E_out ! sampled outgoing energy
real(8), intent(out) :: mu ! sampled scattering cosine
! Sample cosine of scattering angle
if (this%fission) then
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! For fission, the angle is not used, so just assign a dummy value
mu = ONE
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
elseif (allocated(this%angle%energy)) then
mu = this%angle%sample(E_in)
else
! no angle distribution given => assume isotropic for all energies
mu = TWO*prn() - ONE
end if
! Sample outgoing energy
E_out = this%energy%sample(E_in)
end subroutine uncorrelated_sample
end module secondary_uncorrelated

49
src/state_point.F90
View file

@ -310,54 +310,7 @@ contains
call write_dataset(tally_group, "n_score_bins", tally%n_score_bins)
allocate(str_array(size(tally%score_bins)))
do j = 1, size(tally%score_bins)
select case(tally%score_bins(j))
case (SCORE_FLUX)
str_array(j) = "flux"
case (SCORE_TOTAL)
str_array(j) = "total"
case (SCORE_SCATTER)
str_array(j) = "scatter"
case (SCORE_NU_SCATTER)
str_array(j) = "nu-scatter"
case (SCORE_SCATTER_N)
str_array(j) = "scatter-n"
case (SCORE_SCATTER_PN)
str_array(j) = "scatter-pn"
case (SCORE_NU_SCATTER_N)
str_array(j) = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
str_array(j) = "nu-scatter-pn"
case (SCORE_TRANSPORT)
str_array(j) = "transport"
case (SCORE_N_1N)
str_array(j) = "n1n"
case (SCORE_ABSORPTION)
str_array(j) = "absorption"
case (SCORE_FISSION)
str_array(j) = "fission"
case (SCORE_NU_FISSION)
str_array(j) = "nu-fission"
case (SCORE_DELAYED_NU_FISSION)
str_array(j) = "delayed-nu-fission"
case (SCORE_KAPPA_FISSION)
str_array(j) = "kappa-fission"
case (SCORE_CURRENT)
str_array(j) = "current"
case (SCORE_FLUX_YN)
str_array(j) = "flux-yn"
case (SCORE_TOTAL_YN)
str_array(j) = "total-yn"
case (SCORE_SCATTER_YN)
str_array(j) = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
str_array(j) = "nu-scatter-yn"
case (SCORE_EVENTS)
str_array(j) = "events"
case (SCORE_INVERSE_VELOCITY)
str_array(j) = "inverse-velocity"
case default
str_array(j) = reaction_name(tally%score_bins(j))
end select
str_array(j) = reaction_name(tally%score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
call write_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins)

63
src/summary.F90
View file

@ -45,7 +45,7 @@ contains
call write_attribute_string(file_id, "n_batches", &
"description", "Total number of batches")
! Write eigenvalue information
! Write eigenvalue information
if (run_mode == MODE_EIGENVALUE) then
! write number of inactive/active batches and generations/batch
call write_dataset(file_id, "n_inactive", n_inactive)
@ -182,8 +182,10 @@ contains
case (CELL_FILL)
call write_dataset(cell_group, "fill_type", "universe")
call write_dataset(cell_group, "fill", universes(c%fill)%id)
if (size(c%offset) > 0) then
call write_dataset(cell_group, "offset", c%offset)
if (allocated(c%offset)) then
if (size(c%offset) > 0) then
call write_dataset(cell_group, "offset", c%offset)
end if
end if
if (allocated(c%translation)) then
@ -370,8 +372,10 @@ contains
call write_dataset(lattice_group, "outer", lat%outer)
! Write distribcell offsets if present
if (size(lat%offset) > 0) then
call write_dataset(lattice_group, "offsets", lat%offset)
if (allocated(lat%offset)) then
if (size(lat%offset) > 0) then
call write_dataset(lattice_group, "offsets", lat%offset)
end if
end if
select type (lat)
@ -643,54 +647,7 @@ contains
call write_dataset(tally_group, "n_score_bins", t%n_score_bins)
allocate(str_array(size(t%score_bins)))
do j = 1, size(t%score_bins)
select case(t%score_bins(j))
case (SCORE_FLUX)
str_array(j) = "flux"
case (SCORE_TOTAL)
str_array(j) = "total"
case (SCORE_SCATTER)
str_array(j) = "scatter"
case (SCORE_NU_SCATTER)
str_array(j) = "nu-scatter"
case (SCORE_SCATTER_N)
str_array(j) = "scatter-n"
case (SCORE_SCATTER_PN)
str_array(j) = "scatter-pn"
case (SCORE_NU_SCATTER_N)
str_array(j) = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
str_array(j) = "nu-scatter-pn"
case (SCORE_TRANSPORT)
str_array(j) = "transport"
case (SCORE_N_1N)
str_array(j) = "n1n"
case (SCORE_ABSORPTION)
str_array(j) = "absorption"
case (SCORE_FISSION)
str_array(j) = "fission"
case (SCORE_NU_FISSION)
str_array(j) = "nu-fission"
case (SCORE_DELAYED_NU_FISSION)
str_array(j) = "delayed-nu-fission"
case (SCORE_KAPPA_FISSION)
str_array(j) = "kappa-fission"
case (SCORE_CURRENT)
str_array(j) = "current"
case (SCORE_FLUX_YN)
str_array(j) = "flux-yn"
case (SCORE_TOTAL_YN)
str_array(j) = "total-yn"
case (SCORE_SCATTER_YN)
str_array(j) = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
str_array(j) = "nu-scatter-yn"
case (SCORE_EVENTS)
str_array(j) = "events"
case (SCORE_INVERSE_VELOCITY)
str_array(j) = "inverse-velocity"
case default
str_array(j) = reaction_name(t%score_bins(j))
end select
str_array(j) = reaction_name(t%score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)

1
src/tracking.F90
View file

@ -202,6 +202,7 @@ contains
if (p % n_secondary > 0) then
call p % initialize_from_source(p % secondary_bank(p % n_secondary))
p % n_secondary = p % n_secondary - 1
n_event = 0
! Enter new particle in particle track file
if (p % write_track) call add_particle_track()

10
tests/cleanup
View file

@ -1,10 +0,0 @@
#!/bin/bash
# This simple script ensures that all binary
# output files have been deleted in all the
# folders. This can occur if a previous error
# occurred and the test suite was rerun without
# deleting left over binary files. This will
# cause an assertion error in some of the
# tests.
find . \( -name "*.h5" -o -name "*.ppm" \) -exec rm -f {} \;

161
tests/run_tests.py
View file

@ -8,7 +8,7 @@ import shutil
import re
import glob
import socket
from subprocess import call
from subprocess import call, check_output
from collections import OrderedDict
from optparse import OptionParser
@ -42,9 +42,9 @@ parser.add_option("-s", "--script", action="store_true", dest="script",
# Default compiler paths
FC='gfortran'
MPI_DIR='/opt/mpich/3.1.3-gnu'
HDF5_DIR='/opt/hdf5/1.8.15-gnu'
PHDF5_DIR='/opt/phdf5/1.8.15-gnu'
MPI_DIR='/opt/mpich/3.2-gnu'
HDF5_DIR='/opt/hdf5/1.8.16-gnu'
PHDF5_DIR='/opt/phdf5/1.8.16-gnu'
# Script mode for extra capability
script_mode = False
@ -73,11 +73,13 @@ set(CTEST_UPDATE_COMMAND "git")
set(CTEST_CONFIGURE_COMMAND "${{CMAKE_COMMAND}} -H${{CTEST_SOURCE_DIRECTORY}} -B${{CTEST_BINARY_DIRECTORY}} ${{CTEST_BUILD_OPTIONS}}")
set(CTEST_MEMORYCHECK_COMMAND "{valgrind_cmd}")
set(CTEST_MEMORYCHECK_COMMAND_OPTIONS "--tool=memcheck --leak-check=yes --show-reachable=yes --num-callers=20 --track-fds=yes")
set(CTEST_MEMORYCHECK_SUPPRESSIONS_FILE ${{CTEST_SOURCE_DIRECTORY}}/../tests/valgrind.supp)
#set(CTEST_MEMORYCHECK_SUPPRESSIONS_FILE ${{CTEST_SOURCE_DIRECTORY}}/../tests/valgrind.supp)
set(MEM_CHECK {mem_check})
if(MEM_CHECK)
set(ENV{{MEM_CHECK}} ${{MEM_CHECK}})
endif()
set(CTEST_COVERAGE_COMMAND "{gcov_cmd}")
set(CTEST_COVERAGE_COMMAND "gcov")
set(COVERAGE {coverage})
set(ENV{{COVERAGE}} ${{COVERAGE}})
@ -87,9 +89,11 @@ ctest_start("{dashboard}")
ctest_configure(RETURN_VALUE res)
{update}
ctest_build(RETURN_VALUE res)
if(NOT MEM_CHECK)
ctest_test({tests} PARALLEL_LEVEL {n_procs}, RETURN_VALUE res)
endif()
if(MEM_CHECK)
ctest_memcheck({tests}, RETURN_VALUE res)
ctest_memcheck({tests} RETURN_VALUE res)
endif(MEM_CHECK)
if(COVERAGE)
ctest_coverage(RETURN_VALUE res)
@ -105,6 +109,32 @@ endif()
# Define test data structure
tests = OrderedDict()
def cleanup(path):
"""Remove generated output files."""
for dirpath, dirnames, filenames in os.walk(path):
for fname in filenames:
for ext in ['.h5', '.ppm', '.voxel']:
if fname.endswith(ext):
os.remove(os.path.join(dirpath, fname))
def which(program):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
fpath, fname = os.path.split(program)
if fpath:
if is_exe(program):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
exe_file = os.path.join(path, program)
if is_exe(exe_file):
return exe_file
return None
class Test(object):
def __init__(self, name, debug=False, optimize=False, mpi=False, openmp=False,
phdf5=False, valgrind=False, coverage=False):
@ -119,8 +149,6 @@ class Test(object):
self.success = True
self.msg = None
self.skipped = False
self.valgrind_cmd = ""
self.gcov_cmd = ""
self.cmake = ['cmake', '-H..', '-Bbuild',
'-DPYTHON_EXECUTABLE=' + sys.executable]
@ -231,42 +259,6 @@ class Test(object):
self.success = False
self.msg = 'Failed on testing.'
# Checks to see if file exists in PWD or PATH
def check_compiler(self):
result = False
if os.path.isfile(self.fc):
result = True
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, self.fc)):
result = True
if not result:
self.msg = 'Compiler not found: {0}'.\
format((os.path.join(path, self.fc)))
self.success = False
# Get valgrind command from user's environment
def find_valgrind(self):
result = False
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, 'valgrind')):
self.valgrind_cmd = os.path.join(path, 'valgrind')
result = True
break
if not result:
self.msg = 'valgrind not found.'
self.success = False
# Get coverage command from user's environment
def find_coverage(self):
result = False
for path in os.environ["PATH"].split(":"):
if os.path.isfile(os.path.join(path, 'gcov')):
self.gcov_cmd = os.path.join(path, 'gcov')
result = True
break
if not result:
self.msg = 'gcov not found.'
self.success = False
# Simple function to add a test to the global tests dictionary
def add_test(name, debug=False, optimize=False, mpi=False, openmp=False,\
@ -342,7 +334,7 @@ else:
# Setup CTest script vars. Not used in non-script mode
pwd = os.getcwd()
ctest_vars = {
'source_dir': os.path.join(pwd, '..'),
'source_dir': os.path.join(pwd, os.pardir),
'build_dir': os.path.join(pwd, 'build'),
'host_name': socket.gethostname(),
'dashboard': dash,
@ -363,10 +355,10 @@ else:
# Set up default valgrind tests (subset of all tests)
# Currently takes too long to run all the tests with valgrind
# Only used in script mode
valgrind_default_tests = "basic|cmfd_feed|confidence_intervals|\
density_atombcm|eigenvalue_genperbatch|energy_grid|entropy|\
filter_cell|lattice_multiple|output|plot_background|reflective_plane|\
rotation|salphabeta_multiple|score_absorption|seed|source_energy_mono|\
valgrind_default_tests = "cmfd_feed|confidence_intervals|\
density|eigenvalue_genperbatch|energy_grid|entropy|\
lattice_multiple|output|plotreflective_plane|\
rotation|salphabetascore_absorption|seed|source_energy_mono|\
sourcepoint_batch|statepoint_interval|survival_biasing|\
tally_assumesep|translation|uniform_fs|universe|void"
@ -383,7 +375,7 @@ if len(list(tests.keys())) == 0:
# Begin testing
shutil.rmtree('build', ignore_errors=True)
call(['./cleanup']) # removes all binary and hdf5 output files from tests
cleanup('.')
for key in iter(tests):
test = tests[key]
@ -395,29 +387,34 @@ for key in iter(tests):
sys.stdout.flush()
# Verify fortran compiler exists
test.check_compiler()
if not test.success:
if which(test.fc) is None:
self.msg = 'Compiler not found: {0}'.format(test.fc)
self.success = False
continue
# Get valgrind command
# Verify valgrind command exists
if test.valgrind:
test.find_valgrind()
if not test.success:
continue
valgrind_cmd = which('valgrind')
if valgrind_cmd is None:
self.msg = 'No valgrind executable found.'
self.success = False
continue
else:
valgrind_cmd = ''
# Get coverage command
# Verify gcov/lcov exist
if test.coverage:
test.find_coverage()
if not test.success:
continue
if which('gcov') is None:
self.msg = 'No {} executable found.'.format(exe)
self.success = False
continue
# Set test specific CTest script vars. Not used in non-script mode
ctest_vars.update({'build_name' : test.get_build_name()})
ctest_vars.update({'build_opts' : test.get_build_opts()})
ctest_vars.update({'mem_check' : test.valgrind})
ctest_vars.update({'coverage' : test.coverage})
ctest_vars.update({'valgrind_cmd' : test.valgrind_cmd})
ctest_vars.update({'gcov_cmd' : test.gcov_cmd})
ctest_vars.update({'build_name': test.get_build_name()})
ctest_vars.update({'build_opts': test.get_build_opts()})
ctest_vars.update({'mem_check': test.valgrind})
ctest_vars.update({'coverage': test.coverage})
ctest_vars.update({'valgrind_cmd': valgrind_cmd})
# Check for user custom tests
# INCLUDE is a CTest command that allows for a subset
@ -458,7 +455,7 @@ for key in iter(tests):
test.run_ctests()
# Leave build directory
os.chdir('..')
os.chdir(os.pardir)
# Copy over log file
if script_mode:
@ -471,11 +468,37 @@ for key in iter(tests):
logfilename = logfilename + '_{0}.log'.format(test.name)
shutil.copy(logfile[0], logfilename)
# For coverage builds, use lcov to generate HTML output
if test.coverage:
if which('lcov') is None or which('genhtml') is None:
print('No lcov/genhtml command found. '
'Could not generate coverage report.')
else:
shutil.rmtree('coverage', ignore_errors=True)
call(['lcov', '--directory', '.', '--capture',
'--output-file', 'coverage.info'])
call(['genhtml', '--output-directory', 'coverage', 'coverage.info'])
os.remove('coverage.info')
if test.valgrind:
# Copy memcheck output to memcheck directory
shutil.rmtree('memcheck', ignore_errors=True)
os.mkdir('memcheck')
memcheck_out = glob.glob('build/Testing/Temporary/MemoryChecker.*.log')
for fname in memcheck_out:
shutil.copy(fname, 'memcheck/')
# Remove generated XML files
xml_files = check_output(['git', 'ls-files', '.', '--exclude-standard',
'--others']).split()
for f in xml_files:
os.remove(f)
# Clear build directory and remove binary and hdf5 files
shutil.rmtree('build', ignore_errors=True)
if script_mode:
os.remove('ctestscript.run')
call(['./cleanup'])
cleanup('.')
# Print out summary of results
print('\n' + '='*54)

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

@ -0,0 +1 @@
b9b4222c4beea80fe6083590f6b785303d174972d80671fb661bac8e030db6f4a61648240cfad6162799361fc0e08a23c61d31aff844d978528d6dad5b5fbc63

1
tests/test_asymmetric_lattice/results_true.dat Normal file
View file

@ -0,0 +1 @@
b5f96919ca474cd1c9c9d0acde3b8aac4a1cf636443c72a38b6c5a4221a8ce3e90182aaef2f664e44b9175ca257a89db2328b63e19388ee0e5006de4b3d92ce6

129
tests/test_asymmetric_lattice/test_asymmetric_lattice.py Normal file
View file

@ -0,0 +1,129 @@
#!/usr/bin/env python
import os
import sys
import glob
import hashlib
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import openmc
from openmc.source import Source
from openmc.stats import Box
class AsymmetricLatticeTestHarness(PyAPITestHarness):
def _build_inputs(self):
"""Build an axis-asymmetric lattice of fuel assemblies"""
# Build full core geometry from underlying input set
self._input_set.build_default_materials_and_geometry()
# Extract all universes from the full core geometry
geometry = self._input_set.geometry.geometry
all_univs = geometry.get_all_universes()
# Extract universes encapsulating fuel and water assemblies
water = all_univs[7]
fuel = all_univs[8]
# Construct a 3x3 lattice of fuel assemblies
core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202)
core_lat.dimension = (3, 3)
core_lat.lower_left = (-32.13, -32.13)
core_lat.pitch = (21.42, 21.42)
core_lat.universes = [[fuel, water, water],
[fuel, fuel, fuel],
[water, water, water]]
# Create bounding surfaces
min_x = openmc.XPlane(x0=-32.13, boundary_type='reflective')
max_x = openmc.XPlane(x0=+32.13, boundary_type='reflective')
min_y = openmc.YPlane(y0=-32.13, boundary_type='reflective')
max_y = openmc.YPlane(y0=+32.13, boundary_type='reflective')
min_z = openmc.ZPlane(z0=0, boundary_type='reflective')
max_z = openmc.ZPlane(z0=+32.13, boundary_type='reflective')
# Define root universe
root_univ = openmc.Universe(universe_id=0, name='root universe')
root_cell = openmc.Cell(cell_id=1)
root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z
root_cell.fill = core_lat
root_univ.add_cell(root_cell)
# Over-ride geometry in the input set with this 3x3 lattice
self._input_set.geometry.geometry.root_universe = root_univ
# Initialize a "distribcell" filter for the fuel pin cell
distrib_filter = openmc.Filter(type='distribcell', bins=[27])
# Initialize the tallies
tally = openmc.Tally(name='distribcell tally', tally_id=27)
tally.add_filter(distrib_filter)
tally.add_score('nu-fission')
# Initialize the tallies file
tallies_file = openmc.TalliesFile()
tallies_file.add_tally(tally)
# Assign the tallies file to the input set
self._input_set.tallies = tallies_file
# Build default settings
self._input_set.build_default_settings()
# Specify summary output and correct source sampling box
source = Source(space=Box([-32, -32, 0], [32, 32, 32]))
source.space.only_fissionable = True
self._input_set.settings.source = source
self._input_set.settings.output = {'summary': True}
# Write input XML files
self._input_set.export()
def _get_results(self, hash_output=True):
"""Digest info in statepoint and summary and return as a string."""
# Read the statepoint file
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = openmc.StatePoint(statepoint)
# Read the summary file
summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0]
su = openmc.Summary(summary)
sp.link_with_summary(su)
# Extract the tally of interest
tally = sp.get_tally(name='distribcell tally')
# Create a string of all mean, std. dev. values for both tallies
outstr = ''
outstr += ', '.join(map(str, tally.mean.flatten())) + '\n'
outstr += ', '.join(map(str, tally.std_dev.flatten())) + '\n'
# Extract fuel assembly lattices from the summary
all_cells = su.openmc_geometry.get_all_cells()
fuel = all_cells[80].fill
core = all_cells[1].fill
# Append a string of lattice distribcell offsets to the string
outstr += ', '.join(map(str, fuel.offsets.flatten())) + '\n'
outstr += ', '.join(map(str, core.offsets.flatten())) + '\n'
# Hash the results if necessary
if hash_output:
sha512 = hashlib.sha512()
sha512.update(outstr.encode('utf-8'))
outstr = sha512.hexdigest()
return outstr
def _cleanup(self):
super(AsymmetricLatticeTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = AsymmetricLatticeTestHarness('statepoint.10.h5', True)
harness.main()

8
tests/test_basic/geometry.xml
View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/test_basic/materials.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

2
tests/test_basic/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
3.021779E-01 3.813358E-03

14
tests/test_density/geometry.xml Normal file
View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 3"/>
<surface id="2" type="sphere" coeffs="0 0 0 6"/>
<surface id="3" type="sphere" coeffs="0 0 0 9"/>
<surface id="4" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
<cell id="2" material="2" region="1 -2" />
<cell id="3" material="3" region="2 -3" />
<cell id="4" material="4" region="3 -4" />
</geometry>

26
tests/test_density/materials.xml Normal file
View file

@ -0,0 +1,26 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="2">
<density value="4.5e22" units="atom/cm3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="3">
<density value="12.3e3" units="kg/m3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="4">
<density units="sum" />
<nuclide name="U-235" xs="71c" ao="0.3e-2" />
<nuclide name="U-238" xs="71c" ao="0.5e-1" />
<nuclide name="H-1" xs="71c" ao="0.1e-2" />
</material>
</materials>

2
tests/test_density/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.088237E+00 1.999252E-02

0
tests/test_basic/settings.xml → tests/test_density/settings.xml
View file

0
tests/test_basic/test_basic.py → tests/test_density/test_density.py Executable file → Normal file
View file

8
tests/test_density_atombcm/geometry.xml
View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/test_density_atombcm/materials.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

2
tests/test_density_atombcm/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
1.752274E+00 4.032481E-02

11
tests/test_density_atombcm/test_density_atombcm.py
View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

8
tests/test_density_atomcm3/geometry.xml
View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/test_density_atomcm3/materials.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5e22" units="atom/cm3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

2
tests/test_density_atomcm3/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
1.092376E+00 1.759788E-02

16
tests/test_density_atomcm3/settings.xml
View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

11
tests/test_density_atomcm3/test_density_atomcm3.py
View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

8
tests/test_density_kgm3/geometry.xml
View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

9
tests/test_density_kgm3/materials.xml
View file

@ -1,9 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="12.3e3" units="kg/m3" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

2
tests/test_density_kgm3/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
7.994522E-01 1.065745E-02

11
tests/test_density_kgm3/test_density_kgm3.py
View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

8
tests/test_density_sum/geometry.xml
View file

@ -1,8 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!-- Sphere with radius 10 -->
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

11
tests/test_density_sum/materials.xml
View file

@ -1,11 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density units="sum" />
<nuclide name="U-235" xs="71c" ao="0.3e-2" />
<nuclide name="U-238" xs="71c" ao="0.5e-1" />
<nuclide name="H-1" xs="71c" ao="0.1e-2" />
</material>
</materials>

2
tests/test_density_sum/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
3.231215E-01 6.421320E-03

16
tests/test_density_sum/settings.xml
View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

11
tests/test_density_sum/test_density_sum.py
View file

@ -1,11 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

5
tests/test_energy_laws/geometry.xml Normal file
View file

@ -0,0 +1,5 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0 0 0 100" boundary="vacuum"/>
<cell id="1" material="1" region="-1" />
</geometry>

11
tests/test_energy_laws/materials.xml Normal file
View file

@ -0,0 +1,11 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<material id="1">
<density value="20" units="g/cc" />
<nuclide name="U-233" ao="1.0" />
<nuclide name="H-2" ao="1.0" />
<nuclide name="Na-23" ao="1.0" />
<nuclide name="Ta-181" ao="1.0" />
</material>
</materials>

2
tests/test_energy_laws/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
2.130076E+00 1.938907E-03

7
tests/test_density_atombcm/settings.xml → tests/test_energy_laws/settings.xml
View file

@ -1,16 +1,11 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
<space type="point" parameters="0. 0. 0." />
</source>
</settings>

30
tests/test_energy_laws/test_energy_laws.py Normal file
View file

@ -0,0 +1,30 @@
#!/usr/bin/env python
"""The purpose of this test is to provide coverage of energy distributions that
are not covered in other tests. It has a single material with the following
nuclides:
U-233: Only nuclide that has a Watt fission spectrum
H-2: Only nuclide that has an N-body phase space distribution, in this case for
(n,2n)
Na-23: Has an evaporation spectrum and also has reactions that have multiple
angle-energy distributions, so it provides coverage for both of those
situations.
Ta-181: One of a few nuclides that has reactions with Kalbach-Mann distributions
that use linear-linear interpolation.
"""
import glob
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
if __name__ == '__main__':
harness = TestHarness('statepoint.10.*')
harness.main()

1
tests/test_filter_azimuthal/inputs_true.dat
View file

@ -1 +0,0 @@
57d6fd9cb5180c38efd2729a5dea0708cbd5fd0bf7dcf0c9d5c9cef5d818aeab5a926d03e70dedcf1b60d5740938fb3ba80e6ccdb09c661d159c0893da3bd593

76
tests/test_filter_azimuthal/results_true.dat
View file

@ -1,76 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
4.215917E+01
3.561920E+02
4.174788E+01
3.505184E+02
4.603223E+01
4.242918E+02
4.496760E+01
4.075599E+02
4.088099E+01
3.376516E+02
tally 2:
4.157239E+01
3.482158E+02
4.227810E+01
3.613293E+02
4.376107E+01
3.835007E+02
4.644205E+01
4.327195E+02
4.191554E+01
3.522147E+02
tally 3:
4.215917E+01
3.561920E+02
4.174788E+01
3.505184E+02
4.603223E+01
4.242918E+02
4.496402E+01
4.075053E+02
4.088458E+01
3.377000E+02
tally 4:
1.531988E+01
4.816326E+01
9.274393E+00
1.821174E+01
1.595868E+01
5.124238E+01
1.299895E+00
6.417145E-01
1.510024E+01
4.604170E+01
8.533361E+00
1.462765E+01
1.658141E+01
5.595629E+01
1.427417E+00
6.621807E-01
1.683102E+01
5.741400E+01
9.845257E+00
2.028406E+01
1.773179E+01
6.477077E+01
1.536972E+00
6.111079E-01
1.586070E+01
5.360975E+01
9.928220E+00
2.089005E+01
1.737609E+01
6.161847E+01
1.700608E+00
8.439708E-01
1.607027E+01
5.490113E+01
7.569336E+00
1.280955E+01
1.606086E+01
5.308665E+01
9.898901E-01
3.143027E-01

59
tests/test_filter_azimuthal/test_filter_azimuthal.py
View file

@ -1,59 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterAzimuthalTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='azimuthal',
bins=(-3.1416, -1.8850, -0.6283, 0.6283, 1.8850,
3.1416))
tally1 = openmc.Tally(tally_id=1)
tally1.add_filter(filt1)
tally1.add_score('flux')
tally1.estimator = 'tracklength'
tally2 = openmc.Tally(tally_id=2)
tally2.add_filter(filt1)
tally2.add_score('flux')
tally2.estimator = 'analog'
filt3 = openmc.Filter(type='azimuthal', bins=(5,))
tally3 = openmc.Tally(tally_id=3)
tally3.add_filter(filt3)
tally3.add_score('flux')
tally3.estimator = 'tracklength'
mesh = openmc.Mesh(mesh_id=1)
mesh.lower_left = [-182.07, -182.07]
mesh.upper_right = [182.07, 182.07]
mesh.dimension = [2, 2]
filt_mesh = openmc.Filter(type='mesh', bins=(1,))
tally4 = openmc.Tally(tally_id=4)
tally4.add_filter(filt3)
tally4.add_filter(filt_mesh)
tally4.add_score('flux')
tally4.estimator = 'tracklength'
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally1)
self._input_set.tallies.add_tally(tally2)
self._input_set.tallies.add_tally(tally3)
self._input_set.tallies.add_tally(tally4)
self._input_set.tallies.add_mesh(mesh)
super(FilterAzimuthalTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterAzimuthalTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterAzimuthalTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_cell/inputs_true.dat
View file

@ -1 +0,0 @@
f8359184c02fbab5dca5368689a84924066ab1fb09cae575588ceddd696d5461db577498df9959365d89fe933e9b338390e44e362c603c6f2aa5bcf4acc14b20

11
tests/test_filter_cell/results_true.dat
View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
1.767552E+01
6.295417E+01
3.863588E+00
3.013300E+00
5.356594E+01
5.839391E+02

29
tests/test_filter_cell/test_filter_cell.py
View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cell', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterCellTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_cellborn/inputs_true.dat
View file

@ -1 +0,0 @@
8ae662f8881ce8cdec550069c6233c2c91e9a10f7200af6892cf6f2d77712ccfa17895dbd2eee02e6daf3d665c6ed84b29e17d89ff519e70c37b36d75a431d53

11
tests/test_filter_cellborn/results_true.dat
View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
0.000000E+00
0.000000E+00
8.921179E+01
1.601939E+03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

29
tests/test_filter_cellborn/test_filter_cellborn.py
View file

@ -1,29 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterCellbornTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='cellborn', bins=(10, 21, 22, 23))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterCellbornTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterCellbornTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterCellbornTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_delayedgroup/inputs_true.dat
View file

@ -1 +0,0 @@
a7c8ce7ffbc3a7b965d8a3077a4d9132130561afef19047b279b2d23198e248b09664856a092a32394894e19fef7708cebad99b3839d735c4e98ae0c9af58cb7

15
tests/test_filter_delayedgroup/results_true.dat
View file

@ -1,15 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
8.141852E-04
1.337187E-07
4.849156E-03
4.744020E-06
4.460252E-03
4.015453E-06
1.028479E-02
2.136252E-05
5.002274E-03
5.056965E-06
1.974747E-03
7.882970E-07

30
tests/test_filter_delayedgroup/test_filter_delayedgroup.py
View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterDelayedgroupTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='delayedgroup',
bins=(1, 2, 3, 4, 5, 6))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('delayed-nu-fission')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterDelayedgroupTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterDelayedgroupTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterDelayedgroupTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_energy/inputs_true.dat
View file

@ -1 +0,0 @@
51d3e2c43f36712a7b26c5fa26e0e2ca6fb9af205af04f0f8cd44c6b100e36382417c2c63d711e4677ce3c1958d15072727d5fd32424a3f6eb08d1f3b1c7db5a

11
tests/test_filter_energy/results_true.dat
View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.844008E+01
1.619630E+02
4.425619E+01
3.938244E+02
5.527425E+01
6.120383E+02
9.799897E+00
1.957877E+01

30
tests/test_filter_energy/test_filter_energy.py
View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterEnergyTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='energy',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('total')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterEnergyTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterEnergyTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterEnergyTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_energyout/inputs_true.dat
View file

@ -1 +0,0 @@
f0810606c5f947a9fe03bcfc87de3883ce46f59d8603e02ed30f853ebf301b2dc6bdcd109889801ada9e6e0b7be4932efeca97d4beea875af8c8e3ecb7511444

11
tests/test_filter_energyout/results_true.dat
View file

@ -1,11 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.842000E+01
1.620214E+02
4.361000E+01
3.810139E+02
5.297000E+01
5.616595E+02
6.530000E+00
8.828900E+00

30
tests/test_filter_energyout/test_filter_energyout.py
View file

@ -1,30 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterEnergyoutTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt = openmc.Filter(type='energyout',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt)
tally.add_score('scatter')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterEnergyoutTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterEnergyoutTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterEnergyoutTestHarness('statepoint.10.*', True)
harness.main()

1
tests/test_filter_group_transfer/inputs_true.dat
View file

@ -1 +0,0 @@
c4d4334d44956d6dc9abe854a5e9403d7f8a87ffb04a15a3d128e8d18eb4111f46ca277b751e1b0e836d69527502f9abba115a4b2fc64c38da63a9d57968d860

67
tests/test_filter_group_transfer/results_true.dat
View file

@ -1,67 +0,0 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.576000E+01
1.331666E+02
0.000000E+00
0.000000E+00
7.000000E-02
1.300000E-03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.050675E+00
2.274991E-01
0.000000E+00
0.000000E+00
2.070821E+00
8.886068E-01
2.660000E+00
1.422000E+00
0.000000E+00
0.000000E+00
3.897000E+01
3.042635E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.352932E-01
4.705717E-02
0.000000E+00
0.000000E+00
1.018668E+00
2.090017E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.570000E+00
4.182700E+00
0.000000E+00
0.000000E+00
4.968000E+01
4.940534E+02
6.537406E-02
1.230788E-03
0.000000E+00
0.000000E+00
8.678070E-02
2.482037E-03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.290000E+00
2.178900E+00
1.610879E-01
5.883677E-03
6.530000E+00
8.828900E+00
3.151783E-01
2.052521E-02

34
tests/test_filter_group_transfer/test_filter_group_transfer.py
View file

@ -1,34 +0,0 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class FilterGroupTransferTestHarness(PyAPITestHarness):
def _build_inputs(self):
filt1 = openmc.Filter(type='energy',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
filt2 = openmc.Filter(type='energyout',
bins=(0.0, 0.253e-6, 1.0e-3, 1.0, 20.0))
tally = openmc.Tally(tally_id=1)
tally.add_filter(filt1)
tally.add_filter(filt2)
tally.add_score('scatter')
tally.add_score('nu-fission')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_tally(tally)
super(FilterGroupTransferTestHarness, self)._build_inputs()
def _cleanup(self):
super(FilterGroupTransferTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = FilterGroupTransferTestHarness('statepoint.10.*', True)
harness.main()

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