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

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

Browse source
This commit is contained in:
Sterling Harper 2015-12-30 16:37:53 -05:00
commit 295e0034b8
71 changed files with 6954 additions and 3926 deletions
Download patch file Download diff file Expand all files Collapse all files

4
.gitignore vendored
View file

@ -71,4 +71,6 @@ docs/source/pythonapi/examples/*.xml
docs/source/pythonapi/examples/*.png
docs/source/pythonapi/examples/*.xls
docs/source/pythonapi/examples/mgxs
docs/source/pythonapi/examples/tracks
docs/source/pythonapi/examples/tracks
docs/source/pythonapi/examples/fission-rates
docs/source/pythonapi/examples/plots

25
CMakeLists.txt
View file

@ -235,6 +235,14 @@ if(NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/src/xml/fox/.git)
endif()
add_subdirectory(src/xml/fox)
#===============================================================================
# RPATH information
#===============================================================================
# add the automatically determined parts of the RPATH
# which point to directories outside the build tree to the install RPATH
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
#===============================================================================
# Build OpenMC executable
#===============================================================================
@ -278,14 +286,21 @@ target_link_libraries(${program} ${ldflags} ${HDF5_LIBRARIES} fox_dom)
install(TARGETS ${program} RUNTIME DESTINATION bin)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}/copyright")
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)
find_package(PythonInterp)
if(PYTHONINTERP_FOUND)
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} setup.py install
--prefix=${CMAKE_INSTALL_PREFIX}
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})")
if(debian)
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} setup.py install
--root=debian/openmc --install-layout=deb
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})")
else()
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} setup.py install
--prefix=${CMAKE_INSTALL_PREFIX}
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})")
endif()
endif()
#===============================================================================

4
docs/source/conf.py
View file

@ -55,7 +55,7 @@ copyright = u'2011-2015, Massachusetts Institute of Technology'
# The short X.Y version.
version = "0.7"
# The full version, including alpha/beta/rc tags.
release = "0.7.0"
release = "0.7.1"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
@ -200,7 +200,7 @@ latex_elements = {
\usepackage{enumitem}
\usepackage{amsfonts}
\usepackage{amsmath}
\setlistdepth{9}
\setlistdepth{99}
\usepackage{tikz}
\usetikzlibrary{shapes,snakes,shadows,arrows,calc,decorations.markings,patterns,fit,matrix,spy}
\usepackage{fixltx2e}

63
docs/source/publications.rst
View file

@ -26,6 +26,10 @@ Overviews
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>`_
- Daniel J. Kelly, Brian N. Aviles, Paul K. Romano, Bryan R. Herman,
Nicholas E. Horelik, and Benoit Forget, "Analysis of select BEAVRS PWR
benchmark cycle 1 results using MC21 and OpenMC," *Proc. PHYSOR*, Kyoto,
@ -57,13 +61,8 @@ Coupling and Multi-physics
- Bryan R. Herman, Benoit Forget, and Kord Smith, "Progress toward Monte
Carlo-thermal hydraulic coupling using low-order nonlinear diffusion
acceleration methods." In press, *Ann. Nucl. Energy*,
(2014). `<http://dx.doi.org/10.1016/j.anucene/2014.10.029>`_
- Adam G. Nelson and William R. Martin, "Improved Convergence of Monte Carlo
Generated Multi-Group Scattering Moments," *Proc. Int. Conf. Mathematics and
Computational Methods Applied to Nuclear Science and Engineering*, Sun Valley,
Idaho, May 5--9 (2013).
acceleration methods." *Ann. Nucl. Energy*, **84**, 63-72
(2015). `<http://dx.doi.org/10.1016/j.anucene.2014.10.029>`_
- Bryan R. Herman, Benoit Forget, and Kord Smith, "Utilizing CMFD in OpenMC to
Estimate Dominance Ratio and Adjoint," *Trans. Am. Nucl. Soc.*, **109**,
@ -81,19 +80,65 @@ Geometry
Miscellaneous
-------------
- William Boyd, Sterling Harper, and Paul K. Romano, "Equipping OpenMC for the
big data era," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016.
- Qicang Shen, William Boyd, Benoit Forget, and Kord Smith, "Tally precision
triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**,
637-640 (2015).
- Timothy P. Burke, Brian C. Kiedrowski, and William R. Martin, "Flux and
Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*,
**109**, 683-686 (2013).
------------------------------------
Multi-group Cross Section Generation
------------------------------------
- Adam G. Nelson and William R. Martin, "Improved Monte Carlo tallying of
multi-group scattering moments using the NDPP code," *Trans. Am. Nucl. Soc.*,
**113**, 645-648 (2015)
- Adam G. Nelson and William R. Martin, "Improved Monte Carlo tallying of
multi-group scattering moment matrices," *Trans. Am. Nucl. Soc.*, **110**,
217-220 (2014).
- Adam G. Nelson and William R. Martin, "Improved Convergence of Monte Carlo
Generated Multi-Group Scattering Moments," *Proc. Int. Conf. Mathematics and
Computational Methods Applied to Nuclear Science and Engineering*, Sun Valley,
Idaho, May 5--9 (2013).
------------
Nuclear Data
------------
- Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole
for cross section Doppler broadening," *J. Comput. Phys.*, In Press
(2016). `<http://dx.doi.org/10.1016/jcp.2015.08.013>`_
- Colin Josey, Benoit Forget, and Kord Smith, "Windowed multipole sensitivity to
target accuracy of the optimization procedure," *J. Nucl. Sci. Technol.*,
**52**, 987-992 (2015). `<http://dx.doi.org/10.1080/00223131.2015.1035353>`_
- Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith,
"Optimizations of the energy grid search algorithm in continuous-energy Monte
Carlo particle transport codes", *Comput. Phys. Commun.*, **196**, 134-142
(2015). `<http://dx.doi.org/10.1016/j.cpc.2015.05.025>`_
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and
Forrest B. Brown, "Direct, on-the-fly calculation of unresolved resonance
region cross sections in Monte Carlo simulations," *Proc. Joint
Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015).
- Amanda L. Lund, Andrew R. Siegel, Benoit Forget, Colin Josey, and
Paul K. Romano, "Using fractional cascading to accelerate cross section
lookups in Monte Carlo particle transport calculations," *Proc. Joint
Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015).
- Ronald O. Rahaman, Andrew R. Siegel, and Paul K. Romano, "Monte Carlo
performance analysis for varying cross section parameter regimes,"
*Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015).
- Paul K. Romano and Timothy H. Trumbull, "Comparison of algorithms for Doppler
broadening pointwise tabulated cross sections," *Ann. Nucl. Energy*, **75**,
358--364 (2015). `<http://dx.doi.org/10.1016/j.anucene.2014.08.046>`_
@ -114,6 +159,10 @@ Nuclear Data
Parallelism
-----------
- Paul K. Romano, John R. Tramm, and Andrew R. Siegel, "Efficacy of hardware
threading for Monte Carlo particle transport calculations on multi- and
many-core systems," Accepted, *PHYSOR 2016*, Sun Valley, Idaho, May 1-5, 2016.
- David Ozog, Allen D. Malony, and Andrew R. Siegel, "A performance analysis of
SIMD algorithms for Monte Carlo simulations of nuclear reactor cores,"
*Proc. IEEE Int. Parallel and Distributed Processing Symposium*, Hyderabad,

BIN
docs/source/pythonapi/examples/images/mgxs.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 53 KiB

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

File diff suppressed because one or more lines are too long

13
docs/source/pythonapi/examples/mgxs-part-i.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_mgxs_part_i:
=========================
MGXS Part I: Introduction
=========================
.. only:: html
.. notebook:: mgxs-part-i.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

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

File diff suppressed because one or more lines are too long

13
docs/source/pythonapi/examples/mgxs-part-ii.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_mgxs_part_ii:
===============================
MGXS Part II: Advanced Features
===============================
.. only:: html
.. notebook:: mgxs-part-ii.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

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

File diff suppressed because one or more lines are too long

13
docs/source/pythonapi/examples/mgxs-part-iii.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_mgxs_part_iii:
========================
MGXS Part III: Libraries
========================
.. only:: html
.. notebook:: mgxs-part-iii.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

2454
docs/source/pythonapi/examples/multi-group-cross-sections.ipynb
View file

File diff suppressed because one or more lines are too long

11
docs/source/pythonapi/examples/multi-group-cross-sections.rst
View file

@ -1,11 +0,0 @@
====================================
Multi-Group Cross Section Generation
====================================
.. only:: html
.. notebook:: multi-group-cross-sections.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

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

@ -126,7 +126,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's move on to the geometry. This problem will be a square array of fuel pins, which we can use OpenMC's lattice/universe feature for. The basic universe will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces for fuel and clad, as well as the outer bounding surfaces of the problem."
"Now let's move on to the geometry. This problem will be a square array of fuel pins for which we can use OpenMC's lattice/universe feature. The basic universe will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces for fuel and clad, as well as the outer bounding surfaces of the problem."
]
},
{
@ -155,7 +155,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"With the surfaces defined, we can now create cells that are defined by intersections of half-spaces created by the surfaces."
"With the surfaces defined, we can now construct a fuel pin cell from cells that are defined by intersections of half-spaces created by the surfaces."
]
},
{
@ -192,7 +192,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the pin cell universe, we can construct a 17x17 rectangular lattice with a 1.26cm pitch."
"Using the pin cell universe, we can construct a 17x17 rectangular lattice with a 1.26 cm pitch."
]
},
{
@ -242,7 +242,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We now must create a geometry that is assigned a root universe, put the geometry into a geometry file, and export it to XML."
"We now must create a geometry that is assigned a root universe, put the geometry into a `GeometryFile` object, and export it to XML."
]
},
{

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

File diff suppressed because one or more lines are too long

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

@ -363,7 +363,26 @@
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAxQTFRF\n////chIS6YCRTb/E6kGE+wAAAAFiS0dEAIgFHUgAAAAJcEhZcwAAAEgAAABIAEbJaz4AAALKSURB\nVGje7dpLcqQwDAbgHHE2YeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmN\nP+HDhw8fPnz48Kf6VH9G+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4\nzPji99z0/AJ4n1lfvJ6fnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6\npA0wfln+ho/fwgYYn19C/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tN\nDbSGz7T0SBEWw4vLXzbQ6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X5\n8wZaxWd1+fMGiuFvir8bvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV\n873hB8UnM3xzANtf8nb4dwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7\nT/ppARBvp48UwJnelT5SACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4/\n/Jve+fhsH6Ctv7n8PTzjvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V\n32/o9+fl389Xnx+g5x/o+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6\n/4Le/6D3T/D9V67Y/ZsVQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/\ngPs/0P4TtP8F7r9J3AIO9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTu\nf4X7b+H+X7T/+BPuf3aM8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTUtMTAtMjhUMjE6MTU6MDYtMDQ6MDDSDNQdAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE1LTEwLTI4\nVDIxOjE1OjA2LTA0OjAwo1FsoQAAAABJRU5ErkJggg==\n",
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\n",
"AAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\n",
"QYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB98LGQ4UM+6dthcAAALKSURBVGje7dpLcqQwDAbgHHE2\n",
"YeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n",
"+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\n",
"nl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n",
"/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n",
"6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\n",
"vjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\n",
"dwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\n",
"ACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\n",
"vY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n",
"+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\n",
"QBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n",
"9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n",
"8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTUtMTEtMjVUMTQ6MjA6\n",
"NTEtMDg6MDDVsKLDAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE1LTExLTI1VDE0OjIwOjUxLTA4OjAw\n",
"pO0afwAAAABJRU5ErkJggg==\n"
],
"text/plain": [
"<IPython.core.display.Image object>"
]
@ -573,8 +592,8 @@
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.0\n",
" Git SHA1: 21738db07debeabde824c9b955bd3bf0c9a16366\n",
" Date/Time: 2015-10-28 21:15:07\n",
" Git SHA1: 74ffcb447521c968fb64fdaa63e40598783f2fba\n",
" Date/Time: 2015-11-25 14:20:51\n",
" MPI Processes: 1\n",
"\n",
" ===========================================================================\n",
@ -615,13 +634,13 @@
" 11/1 1.07867 1.05536 +/- 0.01277\n",
" 12/1 1.04203 1.05345 +/- 0.01096\n",
" 13/1 1.04482 1.05237 +/- 0.00955\n",
" 14/1 1.04117 1.05113 +/- 0.00852\n",
" 15/1 1.07581 1.05360 +/- 0.00801\n",
" 16/1 1.04235 1.05257 +/- 0.00731\n",
" 17/1 1.02710 1.05045 +/- 0.00701\n",
" 18/1 1.01970 1.04809 +/- 0.00687\n",
" 19/1 1.01022 1.04538 +/- 0.00691\n",
" 20/1 1.01449 1.04332 +/- 0.00675\n",
" 14/1 1.04116 1.05113 +/- 0.00852\n",
" 15/1 1.07569 1.05358 +/- 0.00800\n",
" 16/1 1.04188 1.05252 +/- 0.00732\n",
" 17/1 1.03775 1.05129 +/- 0.00679\n",
" 18/1 0.98462 1.04616 +/- 0.00808\n",
" 19/1 1.08613 1.04902 +/- 0.00801\n",
" 20/1 1.00571 1.04613 +/- 0.00800\n",
" Creating state point statepoint.20.h5...\n",
"\n",
" ===========================================================================\n",
@ -631,27 +650,27 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 6.3800E-01 seconds\n",
" Reading cross sections = 1.3500E-01 seconds\n",
" Total time in simulation = 2.3556E+01 seconds\n",
" Time in transport only = 2.3532E+01 seconds\n",
" Time in inactive batches = 3.1100E+00 seconds\n",
" Time in active batches = 2.0446E+01 seconds\n",
" Total time for initialization = 7.9600E-01 seconds\n",
" Reading cross sections = 2.1200E-01 seconds\n",
" Total time in simulation = 1.8740E+01 seconds\n",
" Time in transport only = 1.8727E+01 seconds\n",
" Time in inactive batches = 2.5970E+00 seconds\n",
" Time in active batches = 1.6143E+01 seconds\n",
" Time synchronizing fission bank = 2.0000E-03 seconds\n",
" Sampling source sites = 1.0000E-03 seconds\n",
" SEND/RECV source sites = 1.0000E-03 seconds\n",
" Time accumulating tallies = 1.0000E-03 seconds\n",
" Total time for finalization = 3.0000E-03 seconds\n",
" Total time elapsed = 2.4210E+01 seconds\n",
" Calculation Rate (inactive) = 4019.29 neutrons/second\n",
" Calculation Rate (active) = 1834.10 neutrons/second\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 2.0000E-03 seconds\n",
" Total time elapsed = 1.9553E+01 seconds\n",
" Calculation Rate (inactive) = 4813.25 neutrons/second\n",
" Calculation Rate (active) = 2322.99 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.03935 +/- 0.00682\n",
" k-effective (Track-length) = 1.04332 +/- 0.00675\n",
" k-effective (Absorption) = 1.03845 +/- 0.00598\n",
" Combined k-effective = 1.04024 +/- 0.00523\n",
" k-effective (Collision) = 1.04597 +/- 0.00663\n",
" k-effective (Track-length) = 1.04613 +/- 0.00800\n",
" k-effective (Absorption) = 1.04087 +/- 0.00627\n",
" Combined k-effective = 1.04322 +/- 0.00570\n",
" Leakage Fraction = 0.00000 +/- 0.00000\n",
"\n"
]
@ -742,7 +761,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -756,10 +775,10 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>total</td>\n",
" <td>(nu-fission / absorption)</td>\n",
" <td>1.040166</td>\n",
" <td>0.009069</td>\n",
" <td> total</td>\n",
" <td> (nu-fission / absorption)</td>\n",
" <td> 1.040687</td>\n",
" <td> 0.010913</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -767,7 +786,7 @@
],
"text/plain": [
" nuclide score mean std. dev.\n",
"0 total (nu-fission / absorption) 1.040166 0.009069"
"0 total (nu-fission / absorption) 1.040687 0.010913"
]
},
"execution_count": 26,
@ -802,7 +821,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -817,19 +836,19 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>total</td>\n",
" <td>absorption</td>\n",
" <td>0.95938</td>\n",
" <td>0.008187</td>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> total</td>\n",
" <td> absorption</td>\n",
" <td> 0.959302</td>\n",
" <td> 0.010033</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total absorption 0.95938 0.008187"
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total absorption 0.959302 0.010033"
]
},
"execution_count": 27,
@ -862,11 +881,12 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>energy [MeV]</th>\n",
" <th>nuclide</th>\n",
" <th>score</th>\n",
" <th>mean</th>\n",
@ -876,18 +896,19 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>total</td>\n",
" <td>nu-fission</td>\n",
" <td>1.090899</td>\n",
" <td>0.010602</td>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> total</td>\n",
" <td> nu-fission</td>\n",
" <td> 1.09103</td>\n",
" <td> 0.012491</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" nuclide score mean std. dev.\n",
"0 total nu-fission 1.090899 0.010602"
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total nu-fission 1.09103 0.012491"
]
},
"execution_count": 28,
@ -921,7 +942,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -937,12 +958,12 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>10000</td>\n",
" <td>total</td>\n",
" <td>absorption</td>\n",
" <td>0.803413</td>\n",
" <td>0.007031</td>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> absorption</td>\n",
" <td> 0.803182</td>\n",
" <td> 0.008664</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -950,7 +971,7 @@
],
"text/plain": [
" energy [MeV] cell nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) 10000 total absorption 0.803413 0.007031"
"0 (0.0e+00 - 6.2e-01) 10000 total absorption 0.803182 0.008664"
]
},
"execution_count": 29,
@ -982,12 +1003,13 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>energy [MeV]</th>\n",
" <th>cell</th>\n",
" <th>nuclide</th>\n",
" <th>score</th>\n",
" <th>mean</th>\n",
@ -997,19 +1019,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>total</td>\n",
" <td>(nu-fission / absorption)</td>\n",
" <td>1.237053</td>\n",
" <td>0.011765</td>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> (nu-fission / absorption)</td>\n",
" <td> 1.237982</td>\n",
" <td> 0.014179</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" energy [MeV] nuclide score mean std. dev.\n",
"0 (0.0e+00 - 6.2e-01) total (nu-fission / absorption) 1.237053 0.011765"
" energy [MeV] cell nuclide score mean \\\n",
"0 (0.0e+00 - 6.2e-01) 10000 total (nu-fission / absorption) 1.237982 \n",
"\n",
" std. dev. \n",
"0 0.014179 "
]
},
"execution_count": 30,
@ -1040,12 +1066,13 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>energy [MeV]</th>\n",
" <th>cell</th>\n",
" <th>nuclide</th>\n",
" <th>score</th>\n",
" <th>mean</th>\n",
@ -1055,22 +1082,23 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>(0.0e+00 - 6.2e-01)</td>\n",
" <td>total</td>\n",
" <td>(((absorption * nu-fission) * absorption) * (n...</td>\n",
" <td>1.040166</td>\n",
" <td>0.019018</td>\n",
" <td> (0.0e+00 - 6.2e-01)</td>\n",
" <td> 10000</td>\n",
" <td> total</td>\n",
" <td> (((absorption * nu-fission) * absorption) * (n...</td>\n",
" <td> 1.040687</td>\n",
" <td> 0.022989</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" energy [MeV] nuclide \\\n",
"0 (0.0e+00 - 6.2e-01) total \n",
" energy [MeV] cell nuclide \\\n",
"0 (0.0e+00 - 6.2e-01) 10000 total \n",
"\n",
" score mean std. dev. \n",
"0 (((absorption * nu-fission) * absorption) * (n... 1.040166 0.019018 "
"0 (((absorption * nu-fission) * absorption) * (n... 1.040687 0.022989 "
]
},
"execution_count": 31,
@ -1118,7 +1146,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1134,100 +1162,100 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>6.657029e-07</td>\n",
" <td>7.377419e-09</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.000001</td>\n",
" <td> 8.078651e-09</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>2.099891e-01</td>\n",
" <td>2.303838e-03</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.209990</td>\n",
" <td> 2.449396e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>3.564204e-01</td>\n",
" <td>3.951669e-03</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.356117</td>\n",
" <td> 4.364366e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>5.555330e-03</td>\n",
" <td>6.101004e-05</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.005555</td>\n",
" <td> 6.495710e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>7.154887e-03</td>\n",
" <td>8.053460e-05</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.007190</td>\n",
" <td> 7.596666e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-238 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>2.277701e-01</td>\n",
" <td>1.079289e-03</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-238 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.227843</td>\n",
" <td> 1.024510e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(nu-fission / flux)</td>\n",
" <td>8.066738e-03</td>\n",
" <td>5.254797e-05</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (nu-fission / flux)</td>\n",
" <td> 0.008086</td>\n",
" <td> 6.251590e-05</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>(U-235 / total)</td>\n",
" <td>(scatter / flux)</td>\n",
" <td>3.366802e-03</td>\n",
" <td>1.647058e-05</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> (U-235 / total)</td>\n",
" <td> (scatter / flux)</td>\n",
" <td> 0.003365</td>\n",
" <td> 1.646663e-05</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" cell energy [MeV] nuclide score \\\n",
"0 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (nu-fission / flux) \n",
"1 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (scatter / flux) \n",
"2 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (nu-fission / flux) \n",
"3 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (scatter / flux) \n",
"4 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (nu-fission / flux) \n",
"5 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (scatter / flux) \n",
"6 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (nu-fission / flux) \n",
"7 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (scatter / flux) \n",
" cell energy [MeV] nuclide score mean \\\n",
"0 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (nu-fission / flux) 0.000001 \n",
"1 10000 (0.0e+00 - 6.3e-07) (U-238 / total) (scatter / flux) 0.209990 \n",
"2 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (nu-fission / flux) 0.356117 \n",
"3 10000 (0.0e+00 - 6.3e-07) (U-235 / total) (scatter / flux) 0.005555 \n",
"4 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (nu-fission / flux) 0.007190 \n",
"5 10000 (6.3e-07 - 2.0e+01) (U-238 / total) (scatter / flux) 0.227843 \n",
"6 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (nu-fission / flux) 0.008086 \n",
"7 10000 (6.3e-07 - 2.0e+01) (U-235 / total) (scatter / flux) 0.003365 \n",
"\n",
" mean std. dev. \n",
"0 6.657029e-07 7.377419e-09 \n",
"1 2.099891e-01 2.303838e-03 \n",
"2 3.564204e-01 3.951669e-03 \n",
"3 5.555330e-03 6.101004e-05 \n",
"4 7.154887e-03 8.053460e-05 \n",
"5 2.277701e-01 1.079289e-03 \n",
"6 8.066738e-03 5.254797e-05 \n",
"7 3.366802e-03 1.647058e-05 "
" std. dev. \n",
"0 8.078651e-09 \n",
"1 2.449396e-03 \n",
"2 4.364366e-03 \n",
"3 6.495710e-05 \n",
"4 7.596666e-05 \n",
"5 1.024510e-03 \n",
"6 6.251590e-05 \n",
"7 1.646663e-05 "
]
},
"execution_count": 33,
@ -1258,11 +1286,11 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 6.65702880e-07]\n",
" [ 3.56420449e-01]]\n",
"[[[ 6.65302296e-07]\n",
" [ 3.56116716e-01]]\n",
"\n",
" [[ 7.15488656e-03]\n",
" [ 8.06673774e-03]]]\n"
" [[ 7.19004460e-03]\n",
" [ 8.08598751e-03]]]\n"
]
}
],
@ -1290,9 +1318,9 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.00555533]]\n",
"[[[ 0.00555516]]\n",
"\n",
" [[ 0.0033668 ]]]\n"
" [[ 0.00336498]]]\n"
]
}
],
@ -1314,8 +1342,8 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[[ 0.22777006]\n",
" [ 0.0033668 ]]]\n"
"[[[ 0.22784316]\n",
" [ 0.00336498]]]\n"
]
}
],
@ -1344,7 +1372,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1360,39 +1388,39 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>U-238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.000002</td>\n",
" <td>1.283958e-08</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> U-238</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.000002</td>\n",
" <td> 1.450189e-08</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>10000</td>\n",
" <td>(0.0e+00 - 6.3e-07)</td>\n",
" <td>U-235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.868553</td>\n",
" <td>6.880390e-03</td>\n",
" <td> 10000</td>\n",
" <td> (0.0e+00 - 6.3e-07)</td>\n",
" <td> U-235</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.870882</td>\n",
" <td> 7.895515e-03</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>U-238</td>\n",
" <td>nu-fission</td>\n",
" <td>0.082149</td>\n",
" <td>8.837250e-04</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> U-238</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.082484</td>\n",
" <td> 8.253437e-04</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>10000</td>\n",
" <td>(6.3e-07 - 2.0e+01)</td>\n",
" <td>U-235</td>\n",
" <td>nu-fission</td>\n",
" <td>0.092618</td>\n",
" <td>5.195308e-04</td>\n",
" <td> 10000</td>\n",
" <td> (6.3e-07 - 2.0e+01)</td>\n",
" <td> U-235</td>\n",
" <td> nu-fission</td>\n",
" <td> 0.092762</td>\n",
" <td> 6.444580e-04</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1400,10 +1428,10 @@
],
"text/plain": [
" cell energy [MeV] nuclide score mean std. dev.\n",
"0 10000 (0.0e+00 - 6.3e-07) U-238 nu-fission 0.000002 1.283958e-08\n",
"1 10000 (0.0e+00 - 6.3e-07) U-235 nu-fission 0.868553 6.880390e-03\n",
"2 10000 (6.3e-07 - 2.0e+01) U-238 nu-fission 0.082149 8.837250e-04\n",
"3 10000 (6.3e-07 - 2.0e+01) U-235 nu-fission 0.092618 5.195308e-04"
"0 10000 (0.0e+00 - 6.3e-07) U-238 nu-fission 0.000002 1.450189e-08\n",
"1 10000 (0.0e+00 - 6.3e-07) U-235 nu-fission 0.870882 7.895515e-03\n",
"2 10000 (6.3e-07 - 2.0e+01) U-238 nu-fission 0.082484 8.253437e-04\n",
"3 10000 (6.3e-07 - 2.0e+01) U-235 nu-fission 0.092762 6.444580e-04"
]
},
"execution_count": 37,
@ -1427,7 +1455,7 @@
{
"data": {
"text/html": [
"<div>\n",
"<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
@ -1443,84 +1471,84 @@
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>10002</td>\n",
" <td>(1.0e-08 - 1.1e-07)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>4.619398</td>\n",
" <td>0.040124</td>\n",
" <td> 10002</td>\n",
" <td> (1.0e-08 - 1.1e-07)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 4.630154</td>\n",
" <td> 0.044512</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>10002</td>\n",
" <td>(1.1e-07 - 1.2e-06)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.030757</td>\n",
" <td>0.011239</td>\n",
" <td> 10002</td>\n",
" <td> (1.1e-07 - 1.2e-06)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.042984</td>\n",
" <td> 0.011429</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>10002</td>\n",
" <td>(1.2e-06 - 1.3e-05)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>1.658488</td>\n",
" <td>0.009777</td>\n",
" <td> 10002</td>\n",
" <td> (1.2e-06 - 1.3e-05)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 1.657517</td>\n",
" <td> 0.008617</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>10002</td>\n",
" <td>(1.3e-05 - 1.4e-04)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>1.853002</td>\n",
" <td>0.007378</td>\n",
" <td> 10002</td>\n",
" <td> (1.3e-05 - 1.4e-04)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 1.863326</td>\n",
" <td> 0.008848</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>10002</td>\n",
" <td>(1.4e-04 - 1.5e-03)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.050773</td>\n",
" <td>0.012484</td>\n",
" <td> 10002</td>\n",
" <td> (1.4e-04 - 1.5e-03)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.043916</td>\n",
" <td> 0.014195</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>10002</td>\n",
" <td>(1.5e-03 - 1.6e-02)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.131759</td>\n",
" <td>0.007821</td>\n",
" <td> 10002</td>\n",
" <td> (1.5e-03 - 1.6e-02)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.134458</td>\n",
" <td> 0.007561</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>10002</td>\n",
" <td>(1.6e-02 - 1.7e-01)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.213710</td>\n",
" <td>0.015159</td>\n",
" <td> 10002</td>\n",
" <td> (1.6e-02 - 1.7e-01)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.209947</td>\n",
" <td> 0.013848</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>10002</td>\n",
" <td>(1.7e-01 - 1.9e+00)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>2.011925</td>\n",
" <td>0.009406</td>\n",
" <td> 10002</td>\n",
" <td> (1.7e-01 - 1.9e+00)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 2.006967</td>\n",
" <td> 0.009368</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>10002</td>\n",
" <td>(1.9e+00 - 2.0e+01)</td>\n",
" <td>H-1</td>\n",
" <td>scatter</td>\n",
" <td>0.371280</td>\n",
" <td>0.003949</td>\n",
" <td> 10002</td>\n",
" <td> (1.9e+00 - 2.0e+01)</td>\n",
" <td> H-1</td>\n",
" <td> scatter</td>\n",
" <td> 0.373895</td>\n",
" <td> 0.002964</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
@ -1528,15 +1556,15 @@
],
"text/plain": [
" cell energy [MeV] nuclide score mean std. dev.\n",
"0 10002 (1.0e-08 - 1.1e-07) H-1 scatter 4.619398 0.040124\n",
"1 10002 (1.1e-07 - 1.2e-06) H-1 scatter 2.030757 0.011239\n",
"2 10002 (1.2e-06 - 1.3e-05) H-1 scatter 1.658488 0.009777\n",
"3 10002 (1.3e-05 - 1.4e-04) H-1 scatter 1.853002 0.007378\n",
"4 10002 (1.4e-04 - 1.5e-03) H-1 scatter 2.050773 0.012484\n",
"5 10002 (1.5e-03 - 1.6e-02) H-1 scatter 2.131759 0.007821\n",
"6 10002 (1.6e-02 - 1.7e-01) H-1 scatter 2.213710 0.015159\n",
"7 10002 (1.7e-01 - 1.9e+00) H-1 scatter 2.011925 0.009406\n",
"8 10002 (1.9e+00 - 2.0e+01) H-1 scatter 0.371280 0.003949"
"0 10002 (1.0e-08 - 1.1e-07) H-1 scatter 4.630154 0.044512\n",
"1 10002 (1.1e-07 - 1.2e-06) H-1 scatter 2.042984 0.011429\n",
"2 10002 (1.2e-06 - 1.3e-05) H-1 scatter 1.657517 0.008617\n",
"3 10002 (1.3e-05 - 1.4e-04) H-1 scatter 1.863326 0.008848\n",
"4 10002 (1.4e-04 - 1.5e-03) H-1 scatter 2.043916 0.014195\n",
"5 10002 (1.5e-03 - 1.6e-02) H-1 scatter 2.134458 0.007561\n",
"6 10002 (1.6e-02 - 1.7e-01) H-1 scatter 2.209947 0.013848\n",
"7 10002 (1.7e-01 - 1.9e+00) H-1 scatter 2.006967 0.009368\n",
"8 10002 (1.9e+00 - 2.0e+01) H-1 scatter 0.373895 0.002964"
]
},
"execution_count": 38,
@ -1569,7 +1597,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
"version": "2.7.10"
}
},
"nbformat": 4,

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

@ -74,7 +74,9 @@ on a given module or class.
examples/post-processing
examples/pandas-dataframes
examples/tally-arithmetic
examples/multi-group-cross-sections
examples/mgxs-part-i
examples/mgxs-part-ii
examples/mgxs-part-iii
.. _Jupyter: https://jupyter.org/
.. _NumPy: http://www.numpy.org/

84
docs/source/releasenotes.rst
View file

@ -1,9 +1,30 @@
.. _releasenotes:
==============================
Release Notes for OpenMC 0.7.0
Release Notes for OpenMC 0.7.1
==============================
This release of OpenMC provides some substantial improvements over version
0.7.0. Non-simple cell regions can now be defined through the ``|`` (union) and
``~`` (complement) operators. Similar changes in the Python API also allow
complex cell regions to be defined. A true secondary particle bank now exists;
this is crucial for photon transport (to be added in the next minor release). A
rich API for multi-group cross section generation has been added via the
``openmc.mgxs`` Python module.
Various improvements to tallies have also been made. It is now possible to
explicitly specify that a collision estimator be used in a tally. A new
``delayedgroup`` filter and ``delayed-nu-fission`` score allow a user to obtain
delayed fission neutron production rates filtered by delayed group. Finally, the
new ``inverse-velocity`` score may be useful for calculating kinetics
parameters.
.. caution:: In previous versions, depending on how OpenMC was compiled binary
output was either given in HDF5 or a flat binary format. With this
version, all binary output is now HDF5 which means you **must**
have HDF5 in order to install OpenMC. Please consult the user's
guide for instructions on how to compile with HDF5.
-------------------
System Requirements
-------------------
@ -17,36 +38,41 @@ the problem at hand (mostly on the number of nuclides in the problem).
New Features
------------
- Complete Python API
- Python 3 compatability for all scripts
- All scripts consistently named openmc-* and installed together
- New 'distribcell' tally filter for repeated cells
- Ability to specify outer lattice universe
- XML input validation utility (openmc-validate-xml)
- Support for hexagonal lattices
- Material union energy grid method
- Tally triggers
- Remove dependence on PETSc
- Significant OpenMP performance improvements
- Support for Fortran 2008 MPI interface
- Use of Travis CI for continuous integration
- Simplifications and improvements to test suite
- Support for complex cell regions (union and complement operators)
- Generic quadric surface type
- Improved handling of secondary particles
- Binary output is now solely HDF5
- ``openmc.mgxs`` Python module enabling multi-group cross section generation
- Collision estimator for tallies
- Delayed fission neutron production tallies with ability to filter by delayed
group
- Inverse velocity tally score
- Performance improvements for binary search
- Performance improvements for reaction rate tallies
---------
Bug Fixes
---------
- b5f712_: Fix bug in spherical harmonics tallies
- e6675b_: Ensure all constants are double precision
- 04e2c1_: Fix potential bug in sample_nuclide routine
- 6121d9_: Fix bugs related to particle track files
- 2f0e89_: Fixes for nuclide specification in tallies
- 299322_: Bug with material filter when void material present
- d74840_: Fix triggers on tallies with multiple filters
- c29a81_: Correctly handle maximum transport energy
- 3edc23_: Fixes in the nu-scatter score
- 629e3b_: Assume unspecified surface coefficients are zero in Python API
- 5dbe8b_: Fix energy filters for openmc-plot-mesh-tally
- ff66f4_: Fixes in the openmc-plot-mesh-tally script
- 441fd4_: Fix bug in kappa-fission score
- 7e5974_: Allow fixed source simulations from Python API
.. _b5f712: https://github.com/mit-crpg/openmc/commit/b5f712
.. _e6675b: https://github.com/mit-crpg/openmc/commit/e6675b
.. _04e2c1: https://github.com/mit-crpg/openmc/commit/04e2c1
.. _6121d9: https://github.com/mit-crpg/openmc/commit/6121d9
.. _2f0e89: https://github.com/mit-crpg/openmc/commit/2f0e89
.. _299322: https://github.com/mit-crpg/openmc/commit/299322
.. _d74840: https://github.com/mit-crpg/openmc/commit/d74840
.. _c29a81: https://github.com/mit-crpg/openmc/commit/c29a81
.. _3edc23: https://github.com/mit-crpg/openmc/commit/3edc23
.. _629e3b: https://github.com/mit-crpg/openmc/commit/629e3b
.. _5dbe8b: https://github.com/mit-crpg/openmc/commit/5dbe8b
.. _ff66f4: https://github.com/mit-crpg/openmc/commit/ff66f4
.. _441fd4: https://github.com/mit-crpg/openmc/commit/441fd4
.. _7e5974: https://github.com/mit-crpg/openmc/commit/7e5974
------------
Contributors
@ -55,13 +81,11 @@ Contributors
This release contains new contributions from the following people:
- `Will Boyd <wbinventor@gmail.com>`_
- `Matt Ellis <mellis13@mit.edu>`_
- `Sterling Harper <sterlingmharper@mit.edu>`_
- `Bryan Herman <bherman@mit.edu>`_
- `Nicholas Horelik <nicholas.horelik@gmail.com>`_
- `Bryan Herman <hermab53@gmail.com>`_
- `Colin Josey <cjosey@mit.edu>`_
- `William Lyu <PaleNeutron@users.noreply.github.com>`_
- `Adam Nelson <nelsonag@umich.edu>`_
- `Paul Romano <paul.k.romano@gmail.com>`_
- `Anthony Scopatz <scopatz@gmail.com>`_
- `Kelly Rowland <kellylynnerowland@gmail.com>`_
- `Sam Shaner <samuelshaner@gmail.com>`_
- `Jon Walsh <walshjon@mit.edu>`_

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

@ -1114,8 +1114,10 @@ Each ``material`` element can have the following attributes or sub-elements:
An element with attributes/sub-elements called ``value`` and ``units``. The
``value`` attribute is the numeric value of the density while the ``units``
can be "g/cm3", "kg/m3", "atom/b-cm", "atom/cm3", or "sum". The "sum" unit
indicates that the density should be calculated as the sum of the atom
fractions for each nuclide in the material. This should not be used in
indicates that values appearing in ``ao`` attributes for ``<nuclide>`` and
``<element>`` sub-elements are to be interpreted as nuclide/element
densities in atom/b-cm, and the total density of the material is taken as
the sum of all nuclides/elements. The "sum" option cannot be used in
conjunction with weight percents.
*Default*: None

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

@ -8,7 +8,7 @@ Installation and Configuration
Installing on Ubuntu with PPA
-----------------------------
For users with Ubuntu 11.10 or later, a binary package for OpenMC is available
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:
@ -28,6 +28,9 @@ Now OpenMC should be recognized within the repository and can be installed:
sudo apt-get install openmc
Binary packages from this PPA may exist for earlier versions of Ubuntu, but they
are no longer supported.
--------------------
Building from Source
--------------------
@ -74,6 +77,12 @@ Prerequisites
You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
.. important::
OpenMC uses various parts of the HDF5 Fortran 2003 API; as such you
must include ``--enable-fortran2003`` or else OpenMC will not be able
to compile.
On Debian derivatives, HDF5 and/or parallel HDF5 can be installed through
the APT package manager:

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

@ -1,3 +1,5 @@
import numpy as np
import openmc
###############################################################################
@ -115,7 +117,7 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('point', [0., 0., 0.])
settings_file.set_source_space('box', np.concatenate(outer_cube.bounding_box))
settings_file.export_to_xml()
###############################################################################

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

@ -1,3 +1,5 @@
import numpy as np
import openmc
###############################################################################
@ -82,5 +84,5 @@ settings_file = openmc.SettingsFile()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.set_source_space('box', np.concatenate(cell.region.bounding_box))
settings_file.export_to_xml()

2
examples/xml/boxes/settings.xml
View file

@ -10,7 +10,7 @@
<!-- Starting source -->
<source>
<space type="point" parameters="0. 0. 0." />
<space type="box" parameters="-10. -10. -10. 10. 10. 10." />
</source>
</settings>

2
openmc/clean_xml.py
View file

@ -1,7 +1,7 @@
def sort_xml_elements(tree):
# Retrieve all children of the root XML node in the tree
elements = tree.getchildren()
elements = list(tree)
# Initialize empty lists for the sorted and comment elements
sorted_elements = []

8
openmc/cross.py
View file

@ -52,7 +52,7 @@ class CrossScore(object):
self.binary_op = binary_op
def __hash__(self):
return hash(str(self))
return hash(repr(self))
def __eq__(self, other):
return str(other) == str(self)
@ -152,7 +152,7 @@ class CrossNuclide(object):
self.binary_op = binary_op
def __hash__(self):
return hash(str(self))
return hash(repr(self))
def __eq__(self, other):
return str(other) == str(self)
@ -309,7 +309,7 @@ class CrossFilter(object):
clone._right_filter = self.right_filter
clone._binary_op = self.binary_op
clone._type = self.type
clone._bins = self.bins
clone._bins = self._bins
clone._num_bins = self.num_bins
clone._stride = self.stride
@ -356,7 +356,7 @@ class CrossFilter(object):
def type(self, filter_type):
if filter_type not in _FILTER_TYPES.values():
msg = 'Unable to set Filter type to "{0}" since it is not one ' \
'of the supported types'.format(type)
'of the supported types'.format(filter_type)
raise ValueError(msg)
self._type = filter_type

2
openmc/element.py
View file

@ -58,7 +58,7 @@ class Element(object):
return not self == other
def __hash__(self):
return hash((self._name, self._xs))
return hash(repr(self))
def __repr__(self):
string = 'Element - {0}\n'.format(self._name)

2
openmc/filter.py
View file

@ -81,7 +81,7 @@ class Filter(object):
return not self == other
def __hash__(self):
return hash((self.type, tuple(self.bins)))
return hash(repr(self))
def __deepcopy__(self, memo):
existing = memo.get(id(self))

46
openmc/material.py
View file

@ -12,17 +12,13 @@ from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.clean_xml import *
# A list of all IDs for all Materials created
MATERIAL_IDS = []
# A static variable for auto-generated Material IDs
AUTO_MATERIAL_ID = 10000
def reset_auto_material_id():
global AUTO_MATERIAL_ID, MATERIAL_IDS
global AUTO_MATERIAL_ID
AUTO_MATERIAL_ID = 10000
MATERIAL_IDS = []
# Units for density supported by OpenMC
@ -83,6 +79,33 @@ class Material(object):
# If specified, this file will be used instead of composition values
self._distrib_otf_file = None
def __eq__(self, other):
if not isinstance(other, Material):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
# FIXME: We cannot compare densities since OpenMC outputs densities
# in atom/b-cm in summary.h5 irregardless of input units, and we
# cannot compute the sum percent in Python since we lack AWR
#elif self.density != other.density:
# return False
#elif self._nuclides != other._nuclides:
# return False
#elif self._elements != other._elements:
# return False
elif self._sab != other._sab:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -165,26 +188,15 @@ class Material(object):
@id.setter
def id(self, material_id):
global AUTO_MATERIAL_ID, MATERIAL_IDS
# If the Material already has an ID, remove it from global list
if hasattr(self, '_id') and self._id is not None:
MATERIAL_IDS.remove(self._id)
if material_id is None:
global AUTO_MATERIAL_ID
self._id = AUTO_MATERIAL_ID
MATERIAL_IDS.append(AUTO_MATERIAL_ID)
AUTO_MATERIAL_ID += 1
else:
check_type('material ID', material_id, Integral)
if material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to "{0}" since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
check_greater_than('material ID', material_id, 0, equality=True)
self._id = material_id
MATERIAL_IDS.append(material_id)
@name.setter
def name(self, name):

3
openmc/mesh.py
View file

@ -189,6 +189,9 @@ class Mesh(object):
cv.check_length('mesh width', width, 2, 3)
self._width = width
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Mesh\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)

18
openmc/mgxs/library.py
View file

@ -67,6 +67,9 @@ class Library(object):
sp_filename : str
The filename of the statepoint with tally data used to the
compute cross sections
keff : Real or None
The combined keff from the statepoint file with tally data used to
compute cross sections (for eigenvalue calculations only)
name : str, optional
Name of the multi-group cross section library. Used as a label to
identify tallies in OpenMC 'tallies.xml' file.
@ -88,6 +91,7 @@ class Library(object):
self._tally_trigger = None
self._all_mgxs = OrderedDict()
self._sp_filename = None
self._keff = None
self.name = name
self.openmc_geometry = openmc_geometry
@ -114,6 +118,7 @@ class Library(object):
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
clone._all_mgxs = self.all_mgxs
clone._sp_filename = self._sp_filename
clone._keff = self._keff
clone._all_mgxs = OrderedDict()
for domain in self.domains:
@ -199,10 +204,15 @@ class Library(object):
def sp_filename(self):
return self._sp_filename
@property
def keff(self):
return self._keff
@openmc_geometry.setter
def openmc_geometry(self, openmc_geometry):
cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
self._openmc_geometry = openmc_geometry
self._opencg_geometry = None
@name.setter
def name(self, name):
@ -361,6 +371,10 @@ class Library(object):
raise ValueError(msg)
self._sp_filename = statepoint._f.filename
self._openmc_geometry = statepoint.summary.openmc_geometry
if statepoint.run_mode == 'k-eigenvalue':
self._keff = statepoint.k_combined[0]
# Load tallies for each MGXS for each domain and mgxs type
for domain in self.domains:
@ -380,9 +394,7 @@ class Library(object):
----------
domain : Material or Cell or Universe or Integral
The material, cell, or universe object of interest (or its ID)
mgxs_type : {'total', 'transport', 'absorption', 'capture', 'fission',
'nu-fission', 'scatter', 'nu-scatter', 'scatter matrix',
'nu-scatter matrix', 'chi'}
mgxs_type : {'total', 'transport', 'absorption', 'capture', 'fission', 'nu-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'chi'}
The type of multi-group cross section object to return
Returns

13
openmc/mgxs/mgxs.py
View file

@ -1225,28 +1225,29 @@ class MGXS(object):
df = df.drop('score', axis=1)
# Override energy groups bounds with indices
groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
groups = np.repeat(groups, self.num_nuclides)
all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
all_groups = np.repeat(all_groups, self.num_nuclides)
if 'energy [MeV]' in df and 'energyout [MeV]' in df:
df.rename(columns={'energy [MeV]': 'group in'}, inplace=True)
in_groups = np.tile(groups, self.num_subdomains)
in_groups = np.tile(all_groups, self.num_subdomains)
in_groups = np.repeat(in_groups, self.num_groups)
df['group in'] = in_groups
df.rename(columns={'energyout [MeV]': 'group out'}, inplace=True)
out_groups = np.tile(groups, self.num_subdomains * self.num_groups)
out_groups = \
np.tile(all_groups, self.num_subdomains * self.num_groups)
df['group out'] = out_groups
columns = ['group in', 'group out']
elif 'energyout [MeV]' in df:
df.rename(columns={'energyout [MeV]': 'group out'}, inplace=True)
in_groups = np.tile(groups, self.num_subdomains)
in_groups = np.tile(all_groups, self.num_subdomains)
df['group out'] = in_groups
columns = ['group out']
elif 'energy [MeV]' in df:
df.rename(columns={'energy [MeV]': 'group in'}, inplace=True)
in_groups = np.tile(groups, self.num_subdomains)
in_groups = np.tile(all_groups, self.num_subdomains)
df['group in'] = in_groups
columns = ['group in']

2
openmc/nuclide.py
View file

@ -61,7 +61,7 @@ class Nuclide(object):
return not self == other
def __hash__(self):
return hash((self._name, self._xs))
return hash(repr(self))
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)

4
openmc/opencg_compatible.py
View file

@ -726,7 +726,7 @@ def get_openmc_cell(opencg_cell):
openmc_cell.fill = get_openmc_material(fill)
if opencg_cell.rotation:
rotation = np.asarray(opencg_cell.rotation, dtype=np.int)
rotation = np.asarray(opencg_cell.rotation, dtype=np.float64)
openmc_cell.rotation = rotation
if opencg_cell.translation:
@ -882,7 +882,7 @@ def get_opencg_lattice(openmc_lattice):
outer = openmc_lattice.outer
if len(pitch) == 2:
new_pitch = np.ones(3, dtype=np.float64)
new_pitch = np.ones(3, dtype=np.float64) * np.inf
new_pitch[:2] = pitch
pitch = new_pitch

54
openmc/region.py
View file

@ -1,6 +1,8 @@
from abc import ABCMeta, abstractmethod
from collections import Iterable
import numpy as np
from openmc.checkvalue import check_type
@ -29,6 +31,17 @@ class Region(object):
def __str__(self):
return ''
def __eq__(self, other):
if not isinstance(other, type(self)):
return False
elif str(self) != str(other):
return False
else:
return True
def __ne__(self, other):
return not self == other
@staticmethod
def from_expression(expression, surfaces):
"""Generate a region given an infix expression.
@ -207,6 +220,8 @@ class Intersection(Region):
----------
nodes : tuple of Region
Regions to take the intersection of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -220,6 +235,16 @@ class Intersection(Region):
def nodes(self):
return self._nodes
@property
def bounding_box(self):
lower_left = np.array([-np.inf, -np.inf, -np.inf])
upper_right = np.array([np.inf, np.inf, np.inf])
for n in self.nodes:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.maximum(lower_left, lower_left_n)
upper_right[:] = np.minimum(upper_right, upper_right_n)
return lower_left, upper_right
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
@ -246,6 +271,8 @@ class Union(Region):
----------
nodes : tuple of Region
Regions to take the union of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -259,6 +286,16 @@ class Union(Region):
def nodes(self):
return self._nodes
@property
def bounding_box(self):
lower_left = np.array([np.inf, np.inf, np.inf])
upper_right = np.array([-np.inf, -np.inf, -np.inf])
for n in self.nodes:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.minimum(lower_left, lower_left_n)
upper_right[:] = np.maximum(upper_right, upper_right_n)
return lower_left, upper_right
@nodes.setter
def nodes(self, nodes):
check_type('nodes', nodes, Iterable, Region)
@ -289,6 +326,8 @@ class Complement(Region):
----------
node : Region
Regions to take the complement of
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -306,3 +345,18 @@ class Complement(Region):
def node(self, node):
check_type('node', node, Region)
self._node = node
@property
def bounding_box(self):
# Use De Morgan's laws to distribute the complement operator so that it
# only applies to surface half-spaces, thus allowing us to calculate the
# bounding box in the usual recursive manner.
if isinstance(self.node, Union):
temp_region = Intersection(*[~n for n in self.node.nodes])
elif isinstance(self.node, Intersection):
temp_region = Union(*[~n for n in self.node.nodes])
elif isinstance(self.node, Complement):
temp_region = self.node.node
else:
temp_region = ~self.node
return temp_region.bounding_box

76
openmc/settings.py
View file

@ -20,6 +20,8 @@ class SettingsFile(object):
Attributes
----------
run_mode : {'eigenvalue' or 'fixed source'}
The type of calculation to perform (default is 'eigenvalue')
batches : int
Number of batches to simulate
generations_per_batch : int
@ -122,7 +124,9 @@ class SettingsFile(object):
"""
def __init__(self):
# Eigenvalue subelement
# Run mode subelement (default is 'eigenvalue')
self._run_mode = 'eigenvalue'
self._batches = None
self._generations_per_batch = None
self._inactive = None
@ -196,9 +200,13 @@ class SettingsFile(object):
self._dd_count_interactions = False
self._settings_file = ET.Element("settings")
self._eigenvalue_subelement = None
self._run_mode_subelement = None
self._source_element = None
@property
def run_mode(self):
return self._run_mode
@property
def batches(self):
return self._batches
@ -399,6 +407,14 @@ class SettingsFile(object):
def dd_count_interactions(self):
return self._dd_count_interactions
@run_mode.setter
def run_mode(self, run_mode):
if 'run_mode' not in ['eigenvalue', 'fixed source']:
msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \
'and "fixed source" are supported."'.format(run_mode)
raise ValueError(msg)
self._run_mode = run_mode
@batches.setter
def batches(self, batches):
check_type('batches', batches, Integral)
@ -861,57 +877,47 @@ class SettingsFile(object):
self._dd_count_interactions = interactions
def _create_eigenvalue_subelement(self):
self._create_particles_subelement()
self._create_batches_subelement()
self._create_inactive_subelement()
self._create_generations_per_batch_subelement()
self._create_keff_trigger_subelement()
def _create_run_mode_subelement(self):
if self.run_mode == 'eigenvalue':
self._run_mode_subelement = \
ET.SubElement(self._settings_file, "eigenvalue")
self._create_particles_subelement()
self._create_batches_subelement()
self._create_inactive_subelement()
self._create_generations_per_batch_subelement()
self._create_keff_trigger_subelement()
else:
if self._run_mode_subelement is None:
self._run_mode_subelement = \
ET.SubElement(self._settings_file, "fixed_source")
self._create_particles_subelement()
self._create_batches_subelement()
def _create_batches_subelement(self):
if self._batches is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "batches")
element = ET.SubElement(self._run_mode_subelement, "batches")
element.text = str(self._batches)
def _create_generations_per_batch_subelement(self):
if self._generations_per_batch is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement,
element = ET.SubElement(self._run_mode_subelement,
"generations_per_batch")
element.text = str(self._generations_per_batch)
def _create_inactive_subelement(self):
if self._inactive is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "inactive")
element = ET.SubElement(self._run_mode_subelement, "inactive")
element.text = str(self._inactive)
def _create_particles_subelement(self):
if self._particles is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "particles")
element = ET.SubElement(self._run_mode_subelement, "particles")
element.text = str(self._particles)
def _create_keff_trigger_subelement(self):
if self._keff_trigger is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "keff_trigger")
element = ET.SubElement(self._run_mode_subelement, "keff_trigger")
for key in self._keff_trigger:
subelement = ET.SubElement(element, key)
@ -1182,10 +1188,10 @@ class SettingsFile(object):
self._settings_file.clear()
self._source_subelement = None
self._trigger_subelement = None
self._eigenvalue_subelement = None
self._run_mode_subelement = None
self._source_element = None
self._create_eigenvalue_subelement()
self._create_run_mode_subelement()
self._create_source_subelement()
self._create_output_subelement()
self._create_statepoint_subelement()

10
openmc/summary.py
View file

@ -567,7 +567,7 @@ class Summary(object):
"""
for index, material in self.materials.items():
if material._id == material_id:
if material.id == material_id:
return material
return None
@ -588,7 +588,7 @@ class Summary(object):
"""
for index, surface in self.surfaces.items():
if surface._id == surface_id:
if surface.id == surface_id:
return surface
return None
@ -609,7 +609,7 @@ class Summary(object):
"""
for index, cell in self.cells.items():
if cell._id == cell_id:
if cell.id == cell_id:
return cell
return None
@ -630,7 +630,7 @@ class Summary(object):
"""
for index, universe in self.universes.items():
if universe._id == universe_id:
if universe.id == universe_id:
return universe
return None
@ -651,7 +651,7 @@ class Summary(object):
"""
for index, lattice in self.lattices.items():
if lattice._id == lattice_id:
if lattice.id == lattice_id:
return lattice
return None

281
openmc/surface.py
View file

@ -3,6 +3,8 @@ from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
import numpy as np
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.region import Region
@ -136,6 +138,33 @@ class Surface(object):
check_value('boundary type', boundary_type, _BC_TYPES)
self._boundary_type = boundary_type
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. If the half-space is
unbounded in a particular direction, numpy.inf is used to represent
infinity.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
def create_xml_subelement(self):
element = ET.Element("surface")
element.set("id", str(self._id))
@ -194,6 +223,10 @@ class Plane(Surface):
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
self._coeffs['A'] = 1.
self._coeffs['B'] = 0.
self._coeffs['C'] = 0.
self._coeffs['D'] = 0.
if A is not None:
self.a = A
@ -276,6 +309,7 @@ class XPlane(Plane):
self._type = 'x-plane'
self._coeff_keys = ['x0']
self._coeffs['x0'] = 0.
if x0 is not None:
self.x0 = x0
@ -289,6 +323,37 @@ class XPlane(Plane):
check_type('x0 coefficient', x0, Real)
self._coeffs['x0'] = x0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the x-plane surface, the
half-spaces are unbounded in their y- and z- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([self.x0, np.inf, np.inf]))
elif side == '+':
return (np.array([self.x0, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class YPlane(Plane):
"""A plane perpendicular to the y axis, i.e. a surface of the form :math:`y -
@ -322,6 +387,7 @@ class YPlane(Plane):
self._type = 'y-plane'
self._coeff_keys = ['y0']
self._coeffs['y0'] = 0.
if y0 is not None:
self.y0 = y0
@ -335,6 +401,37 @@ class YPlane(Plane):
check_type('y0 coefficient', y0, Real)
self._coeffs['y0'] = y0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the y-plane surface, the
half-spaces are unbounded in their x- and z- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, self.y0, np.inf]))
elif side == '+':
return (np.array([-np.inf, self.y0, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class ZPlane(Plane):
"""A plane perpendicular to the z axis, i.e. a surface of the form :math:`z -
@ -368,6 +465,7 @@ class ZPlane(Plane):
self._type = 'z-plane'
self._coeff_keys = ['z0']
self._coeffs['z0'] = 0.
if z0 is not None:
self.z0 = z0
@ -381,6 +479,37 @@ class ZPlane(Plane):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the z-plane surface, the
half-spaces are unbounded in their x- and y- directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, self.z0]))
elif side == '+':
return (np.array([-np.inf, -np.inf, self.z0]),
np.array([np.inf, np.inf, np.inf]))
class Cylinder(Surface):
"""A cylinder whose length is parallel to the x-, y-, or z-axis.
@ -415,6 +544,7 @@ class Cylinder(Surface):
super(Cylinder, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['R']
self._coeffs['R'] = 1.
if R is not None:
self.r = R
@ -468,6 +598,8 @@ class XCylinder(Cylinder):
self._type = 'x-cylinder'
self._coeff_keys = ['y0', 'z0', 'R']
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
if y0 is not None:
self.y0 = y0
@ -493,6 +625,38 @@ class XCylinder(Cylinder):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the x-cylinder surface,
the negative half-space is unbounded in the x- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]),
np.array([np.inf, self.y0 + self.r, self.z0 + self.r]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class YCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the y-axis. This is a
@ -533,6 +697,8 @@ class YCylinder(Cylinder):
self._type = 'y-cylinder'
self._coeff_keys = ['x0', 'z0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['z0'] = 0.
if x0 is not None:
self.x0 = x0
@ -558,6 +724,38 @@ class YCylinder(Cylinder):
check_type('z0 coefficient', z0, Real)
self._coeffs['z0'] = z0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the y-cylinder surface,
the negative half-space is unbounded in the y- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]),
np.array([self.x0 + self.r, np.inf, self.z0 + self.r]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class ZCylinder(Cylinder):
"""An infinite cylinder whose length is parallel to the z-axis. This is a
@ -598,6 +796,8 @@ class ZCylinder(Cylinder):
self._type = 'z-cylinder'
self._coeff_keys = ['x0', 'y0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
if x0 is not None:
self.x0 = x0
@ -623,6 +823,38 @@ class ZCylinder(Cylinder):
check_type('y0 coefficient', y0, Real)
self._coeffs['y0'] = y0
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. For the z-cylinder surface,
the negative half-space is unbounded in the z- direction and the
positive half-space is unbounded in all directions. To represent
infinity, numpy.inf is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]),
np.array([self.x0 + self.r, self.y0 + self.r, np.inf]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class Sphere(Surface):
"""A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`.
@ -667,6 +899,10 @@ class Sphere(Surface):
self._type = 'sphere'
self._coeff_keys = ['x0', 'y0', 'z0', 'R']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
self._coeffs['R'] = 1.
if x0 is not None:
self.x0 = x0
@ -716,6 +952,39 @@ class Sphere(Surface):
check_type('R coefficient', R, Real)
self._coeffs['R'] = R
def bounding_box(self, side):
"""Determine an axis-aligned bounding box.
An axis-aligned bounding box for surface half-spaces is represented by
its lower-left and upper-right coordinates. The positive half-space of a
sphere is unbounded in all directions. To represent infinity, numpy.inf
is used.
Parameters
----------
side : {'+', '-'}
Indicates the negative or positive half-space
Returns
-------
numpy.array
Lower-left coordinates of the axis-aligned bounding box for the
desired half-space
numpy.array
Upper-right coordinates of the axis-aligned bounding box for the
desired half-space
"""
if side == '-':
return (np.array([self.x0 - self.r, self.y0 - self.r,
self.z0 - self.r]),
np.array([self.x0 + self.r, self.y0 + self.r,
self.z0 + self.r]))
elif side == '+':
return (np.array([-np.inf, -np.inf, -np.inf]),
np.array([np.inf, np.inf, np.inf]))
class Cone(Surface):
"""A conical surface parallel to the x-, y-, or z-axis.
@ -761,6 +1030,10 @@ class Cone(Surface):
super(Cone, self).__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
self._coeffs['x0'] = 0.
self._coeffs['y0'] = 0.
self._coeffs['z0'] = 0.
self._coeffs['R2'] = 1.
if x0 is not None:
self.x0 = x0
@ -982,6 +1255,8 @@ class Quadric(Surface):
self._type = 'quadric'
self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']
for key in self._coeff_keys:
self._coeffs[key] = 0.
if a is not None:
self.a = a
@ -1127,6 +1402,8 @@ class Halfspace(Region):
Surface which divides Euclidean space.
side : {'+', '-'}
Indicates whether the positive or negative half-space is used.
bounding_box : tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
"""
@ -1155,6 +1432,10 @@ class Halfspace(Region):
check_value('side', side, ('+', '-'))
self._side = side
@property
def bounding_box(self):
return self.surface.bounding_box(self.side)
def __str__(self):
return '-' + str(self.surface.id) if self.side == '-' \
else str(self.surface.id)

195
openmc/tallies.py
View file

@ -204,24 +204,7 @@ class Tally(object):
return not self == other
def __hash__(self):
hashable = []
for filter in self.filters:
hashable.append((filter.type, tuple(filter.bins)))
for nuclide in self.nuclides:
hashable.append(nuclide.name)
for score in self.scores:
hashable.append(score)
hashable.append(self.estimator)
hashable.append(self.name)
hashable.append(self._diff_variable)
hashable.append(self._diff_material)
hashable.append(self._diff_nuclide)
return hash(tuple(hashable))
return hash(repr(self))
def __repr__(self):
string = 'Tally\n'
@ -1572,19 +1555,31 @@ class Tally(object):
new_name = '({0} {1} {2})'.format(self.name, binary_op, other.name)
new_tally.name = new_name
# Create copies of self and other tallies to rearrange for tally
# arithmetic
self_copy = copy.deepcopy(self)
other_copy = copy.deepcopy(other)
# Find any shared filters between the two tallies
self_filters = set(self.filters)
other_filters = set(other.filters)
filter_intersect = self_filters.intersection(other_filters)
filter_intersect = []
for filter in self_copy.filters:
if filter in other_copy.filters:
filter_intersect.append(filter)
# Align the shared filters to follow in each tally operand
# Align the shared filters in successive order
for i, filter in enumerate(filter_intersect):
self_index = self.filters.index(filter)
other_filter = other.filters[self_index]
if other_filter != filter:
other = other.swap_filters(filter, other_filter)
self_index = self_copy.filters.index(filter)
other_index = other_copy.filters.index(filter)
data = self._align_tally_data(other)
# If necessary, swap self filter
if self_index != i:
self_copy.swap_filters(filter, self_copy.filters[i], inplace=True)
# If necessary, swap other filter
if other_index != i:
other_copy.swap_filters(filter, other_copy.filters[i], inplace=True)
data = self_copy._align_tally_data(other_copy)
if binary_op == '+':
new_tally._mean = data['self']['mean'] + data['other']['mean']
@ -1615,16 +1610,16 @@ class Tally(object):
new_tally._std_dev = np.abs(new_tally.mean) * \
np.sqrt(first_term**2 + second_term**2)
if self.estimator == other.estimator:
new_tally.estimator = self.estimator
if self.with_summary and other.with_summary:
new_tally.with_summary = self.with_summary
if self.num_realizations == other.num_realizations:
new_tally.num_realizations = self.num_realizations
if self_copy.estimator == other_copy.estimator:
new_tally.estimator = self_copy.estimator
if self_copy.with_summary and other_copy.with_summary:
new_tally.with_summary = self_copy.with_summary
if self_copy.num_realizations == other_copy.num_realizations:
new_tally.num_realizations = self_copy.num_realizations
# If filters are identical, simply reuse them in derived tally
if self.filters == other.filters:
for self_filter in self.filters:
if self_copy.filters == other_copy.filters:
for self_filter in self_copy.filters:
new_tally.add_filter(self_filter)
# Generate filter "outer products" for non-identical filters
@ -1632,24 +1627,24 @@ class Tally(object):
# Find the common longest sequence of shared filters
match = 0
for self_filter, other_filter in zip(self.filters, other.filters):
for self_filter, other_filter in zip(self_copy.filters, other_copy.filters):
if self_filter == other_filter:
match += 1
else:
break
match_filters = self.filters[:match]
cross_filters = [self.filters[match:], other.filters[match:]]
match_filters = self_copy.filters[:match]
cross_filters = [self_copy.filters[match:], other_copy.filters[match:]]
# Simply reuse shared filters in derived tally
for filter in match_filters:
new_tally.add_filter(filter)
# Use cross filters to combine non-shared filters in derived tally
if len(self.filters) != match and len(other.filters) == match:
if len(self_copy.filters) != match and len(other_copy.filters) == match:
for filter in cross_filters[0]:
new_tally.add_filter(filter)
elif len(other.filters) == match and len(other.filters) != match:
elif len(self_copy.filters) == match and len(other_copy.filters) != match:
for filter in cross_filters[1]:
new_tally.add_filter(filter)
else:
@ -1658,23 +1653,23 @@ class Tally(object):
new_tally.add_filter(new_filter)
# Generate score "outer products"
if self.scores == other.scores:
new_tally.num_score_bins = self.num_score_bins
for self_score in self.scores:
if self_copy.scores == other_copy.scores:
new_tally.num_score_bins = self_copy.num_score_bins
for self_score in self_copy.scores:
new_tally.add_score(self_score)
else:
new_tally.num_score_bins = self.num_score_bins * other.num_score_bins
all_scores = [self.scores, other.scores]
new_tally.num_score_bins = self_copy.num_score_bins * other_copy.num_score_bins
all_scores = [self_copy.scores, other_copy.scores]
for self_score, other_score in itertools.product(*all_scores):
new_score = CrossScore(self_score, other_score, binary_op)
new_tally.add_score(new_score)
# Generate nuclide "outer products"
if self.nuclides == other.nuclides:
for self_nuclide in self.nuclides:
if self_copy.nuclides == other_copy.nuclides:
for self_nuclide in self_copy.nuclides:
new_tally.nuclides.append(self_nuclide)
else:
all_nuclides = [self.nuclides, other.nuclides]
all_nuclides = [self_copy.nuclides, other_copy.nuclides]
for self_nuclide, other_nuclide in itertools.product(*all_nuclides):
new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op)
new_tally.add_nuclide(new_nuclide)
@ -1734,8 +1729,8 @@ class Tally(object):
self_repeat_factor *= filter.num_bins
# Tile / repeat the tally data for the tally outer product
self_shape = list(self.mean.shape)
other_shape = list(other.mean.shape)
self_shape = list(self_mean.shape)
other_shape = list(other_mean.shape)
self_shape[0] *= self_repeat_factor
self_mean = np.repeat(self_mean, self_repeat_factor)
self_std_dev = np.repeat(self_std_dev, self_repeat_factor)
@ -1743,7 +1738,8 @@ class Tally(object):
if self_repeat_factor == 1:
other_shape[0] *= other_tile_factor
other_mean = np.repeat(other_mean, other_tile_factor, axis=0)
other_std_dev = np.repeat(other_std_dev, other_tile_factor, axis=0)
other_std_dev = np.repeat(other_std_dev, other_tile_factor,
axis=0)
else:
other_mean = np.tile(other_mean, (other_tile_factor, 1, 1))
other_std_dev = np.tile(other_std_dev, (other_tile_factor, 1, 1))
@ -1762,7 +1758,11 @@ class Tally(object):
self_repeat_factor = other.num_nuclides
other_tile_factor = self.num_nuclides
# Replicate the data
# Tile / repeat the tally data for the tally outer product
self_shape = list(self_mean.shape)
other_shape = list(other_mean.shape)
self_shape[1] *= self_repeat_factor
other_shape[1] *= other_tile_factor
self_mean = np.repeat(self_mean, self_repeat_factor, axis=1)
other_mean = np.tile(other_mean, (1, other_tile_factor, 1))
self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=1)
@ -1770,10 +1770,10 @@ class Tally(object):
# NumPy repeat and tile routines return 1D flattened arrays
# Reshape arrays as 3D with filters, nuclides and scores axes
self_shape = list(self.mean.shape)
self_shape[1] *= self_repeat_factor
self_mean.shape = tuple(self_shape)
self_std_dev.shape = tuple(self_shape)
other_mean.shape = tuple(other_shape)
other_std_dev.shape = tuple(other_shape)
if self.scores != other.scores:
@ -1782,7 +1782,11 @@ class Tally(object):
self_repeat_factor = other.num_score_bins
other_tile_factor = self.num_score_bins
# Replicate the data
# Tile / repeat the tally data for the tally outer product
self_shape = list(self_mean.shape)
other_shape = list(other_mean.shape)
self_shape[2] *= self_repeat_factor
other_shape[2] *= other_tile_factor
self_mean = np.repeat(self_mean, self_repeat_factor, axis=2)
other_mean = np.tile(other_mean, (1, 1, other_tile_factor))
self_std_dev = np.repeat(self_std_dev, self_repeat_factor, axis=2)
@ -1790,10 +1794,10 @@ class Tally(object):
# NumPy repeat and tile routines return 1D flattened arrays
# Reshape arrays as 3D with filters, nuclides and scores axes
self_shape = list(self.mean.shape)
self_shape[2] *= self_repeat_factor
self_mean.shape = tuple(self_shape)
self_std_dev.shape = tuple(self_shape)
other_mean.shape = tuple(other_shape)
other_std_dev.shape = tuple(other_shape)
data = {}
data['self'] = {}
@ -1804,7 +1808,7 @@ class Tally(object):
data['other']['std. dev.'] = other_std_dev
return data
def swap_filters(self, filter1, filter2):
def swap_filters(self, filter1, filter2, inplace=False):
"""Reverse the ordering of two filters in this tally
This is a helper method for tally arithmetic which helps align the data
@ -1819,10 +1823,15 @@ class Tally(object):
filter2 : Filter
The filter to swap with filter1
inplace : bool, optional
Whether to perform operation inplace or return new tally with the
filters swapped.
Returns
-------
swap_tally
A copy of this tally with the filters swapped
If inplace is false, a copy of this tally with the filters swapped.
Otherwise, nothing is returned.
Raises
------
@ -1853,7 +1862,15 @@ class Tally(object):
'does not contain such a filter'.format(filter2.type, self.id)
raise ValueError(msg)
swap_tally = copy.deepcopy(self)
# Create a copy of the tally that preserves the original data formatting
# throughout swapping process
tally_copy = copy.deepcopy(self)
# Set the swap tally
if inplace:
swap_tally = self
else:
swap_tally = copy.deepcopy(self)
# Swap the filters in the copied version of this Tally
filter1_index = swap_tally.filters.index(filter1)
@ -1880,42 +1897,43 @@ class Tally(object):
filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)]
# Adjust the sum data array to relect the new filter order
if self.sum is not None:
if swap_tally.sum is not None:
for bin1, bin2 in itertools.product(filter1_bins, filter2_bins):
filter_bins = [(bin1,), (bin2,)]
data = self.get_values(filters=filters,
filter_bins=filter_bins, value='sum')
data = tally_copy.get_values(
filters=filters, filter_bins=filter_bins, value='sum')
indices = swap_tally.get_filter_indices(filters, filter_bins)
swap_tally.sum[indices, :, :] = data
# Adjust the sum_sq data array to relect the new filter order
if self.sum_sq is not None:
if swap_tally.sum_sq is not None:
for bin1, bin2 in itertools.product(filter1_bins, filter2_bins):
filter_bins = [(bin1,), (bin2,)]
data = self.get_values(filters=filters,
filter_bins=filter_bins, value='sum_sq')
data = tally_copy.get_values(
filters=filters, filter_bins=filter_bins, value='sum_sq')
indices = swap_tally.get_filter_indices(filters, filter_bins)
swap_tally.sum_sq[indices, :, :] = data
# Adjust the mean data array to relect the new filter order
if self.mean is not None:
if swap_tally.mean is not None:
for bin1, bin2 in itertools.product(filter1_bins, filter2_bins):
filter_bins = [(bin1,), (bin2,)]
data = self.get_values(filters=filters,
filter_bins=filter_bins, value='mean')
data = tally_copy.get_values(
filters=filters, filter_bins=filter_bins, value='mean')
indices = swap_tally.get_filter_indices(filters, filter_bins)
swap_tally._mean[indices, :, :] = data
# Adjust the std_dev data array to relect the new filter order
if self.std_dev is not None:
if swap_tally.std_dev is not None:
for bin1, bin2 in itertools.product(filter1_bins, filter2_bins):
filter_bins = [(bin1,), (bin2,)]
data = self.get_values(filters=filters,
filter_bins=filter_bins, value='std_dev')
data = tally_copy.get_values(
filters=filters, filter_bins=filter_bins, value='std_dev')
indices = swap_tally.get_filter_indices(filters, filter_bins)
swap_tally._std_dev[indices, :, :] = data
return swap_tally
if not inplace:
return swap_tally
def __add__(self, other):
"""Adds this tally to another tally or scalar value.
@ -2505,7 +2523,7 @@ class Tally(object):
return new_tally
def summation(self, scores=[], filter_type=None,
filter_bins=[], nuclides=[]):
filter_bins=[], nuclides=[], remove_filter=False):
"""Vectorized sum of tally data across scores, filter bins and/or
nuclides using tally addition.
@ -2534,6 +2552,9 @@ class Tally(object):
nuclides : list of str
A list of nuclide name strings to sum across
(e.g., ['U-235', 'U-238']; default is [])
remove_filter : bool
If a filter is being summed over, this bool indicates whether to
remove that filter in the returned tally. Default is False.
Returns
-------
@ -2557,7 +2578,14 @@ class Tally(object):
# Sum across any filter bins specified by the user
if filter_type in _FILTER_TYPES:
filter_bins = [[(filter_bin,)] for filter_bin in filter_bins]
# If user did not specify filter bins, sum across all bins
if len(filter_bins) == 0:
filter = self.find_filter(filter_type)
filter_bins = [[(filter.get_bin(i),)] for i in range(filter.num_bins)]
else:
filter_bins = [[(filter_bin,)] for filter_bin in filter_bins]
filters = [[filter_type]]
# If user did not specify a filter type, do not sum across filter bins
else:
@ -2582,12 +2610,17 @@ class Tally(object):
# Accumulate this Tally slice into the Tally sum
tally_sum += tally_slice
# Add back the filter(s) which were summed across to derived tally
for filter_type in summed_filters:
filters = summed_filters[filter_type]
for i in range(1, len(filters)):
filters[i] = CrossFilter(filters[i-1], filters[i], '+')
tally_sum.add_filter(filters[-1])
# Add back the filter(s) which were summed across to derived tally,
# if filter bins were input; otherwise, leave out summed filter(s)
if remove_filter and filter_type is not None:
# Rename tally sum indicating a summation over a particular filter
tally_sum.name = 'sum({0}, {1})'.format(self.name, filter_type)
else:
for summed_filter_type in summed_filters:
filters = summed_filters[summed_filter_type]
for i in range(1, len(filters)):
filters[i] = CrossFilter(filters[i-1], filters[i], '+')
tally_sum.add_filter(filters[-1])
return tally_sum

128
openmc/universe.py
View file

@ -15,6 +15,10 @@ from openmc.region import Region, Intersection, Complement
if sys.version_info[0] >= 3:
basestring = str
# DeprecationWarning filter for the Cell.add_surface(...) method
warnings.simplefilter('always', DeprecationWarning)
# A static variable for auto-generated Cell IDs
AUTO_CELL_ID = 10000
@ -73,6 +77,30 @@ class Cell(object):
self._translation = None
self._offsets = None
def __eq__(self, other):
if not isinstance(other, Cell):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif self.fill != other.fill:
return False
elif self.region != other.region:
return False
elif self.rotation != other.rotation:
return False
elif self.translation != other.translation:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Cell\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -216,7 +244,6 @@ class Cell(object):
"""
warnings.simplefilter('always', DeprecationWarning)
warnings.warn("Cell.add_surface(...) has been deprecated and may be "
"removed in a future version. The region for a Cell "
"should be defined using the region property directly.",
@ -443,6 +470,34 @@ class Universe(object):
self._cell_offsets = OrderedDict()
self._num_regions = 0
def __eq__(self, other):
if not isinstance(other, Universe):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif self.cells != other.cells:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'Universe\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tCells', '=\t',
list(self._cells.keys()))
string += '{0: <16}{1}{2}\n'.format('\t# Regions', '=\t',
self._num_regions)
return string
@property
def id(self):
return self._id
@ -630,16 +685,6 @@ class Universe(object):
return universes
def __repr__(self):
string = 'Universe\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
string += '{0: <16}{1}{2}\n'.format('\tCells', '=\t',
list(self._cells.keys()))
string += '{0: <16}{1}{2}\n'.format('\t# Regions', '=\t',
self._num_regions)
return string
def create_xml_subelement(self, xml_element):
# Iterate over all Cells
@ -695,6 +740,25 @@ class Lattice(object):
self._outer = None
self._universes = None
def __eq__(self, other):
if not isinstance(other, Lattice):
return False
elif self.id != other.id:
return False
elif self.name != other.name:
return False
elif self.pitch != other.pitch:
return False
elif self.outer != other.outer:
return False
elif self.universes != other.universes:
return False
else:
return True
def __ne__(self, other):
return not self == other
@property
def id(self):
return self._id
@ -894,6 +958,24 @@ class RectLattice(Lattice):
self._lower_left = None
self._offsets = None
def __eq__(self, other):
if not isinstance(other, RectLattice):
return False
elif not super(RectLattice, self).__eq__(other):
return False
elif self.dimension != other.dimension:
return False
elif self.lower_left != other.lower_left:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'RectLattice\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -1037,7 +1119,7 @@ class RectLattice(Lattice):
for z in range(self._dimension[2]):
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[x][y][z]
universe = self._universes[z][y][x]
# Append Universe ID to the Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
@ -1055,7 +1137,7 @@ class RectLattice(Lattice):
else:
for y in range(self._dimension[1]):
for x in range(self._dimension[0]):
universe = self._universes[x][y]
universe = self._universes[y][x]
# Append Universe ID to Lattice XML subelement
universe_ids += '{0} '.format(universe._id)
@ -1111,6 +1193,26 @@ class HexLattice(Lattice):
self._num_axial = None
self._center = None
def __eq__(self, other):
if not isinstance(other, HexLattice):
return False
elif not super(HexLattice, self).__eq__(other):
return False
elif self.num_rings != other.num_rings:
return False
elif self.num_axial != other.num_axial:
return False
elif self.center != other.center:
return False
else:
return True
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash(repr(self))
def __repr__(self):
string = 'HexLattice\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)

2
setup.py
View file

@ -10,7 +10,7 @@ except ImportError:
have_setuptools = False
kwargs = {'name': 'openmc',
'version': '0.7.0',
'version': '0.7.1',
'packages': ['openmc', 'openmc.mgxs'],
'scripts': glob.glob('scripts/openmc-*'),

404
src/ace.F90
View file

@ -45,9 +45,9 @@ contains
integer :: temp_table ! temporary value for sorting
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(Material), pointer :: mat
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(SetChar) :: already_read
! allocate arrays for ACE table storage and cross section cache
@ -234,7 +234,6 @@ contains
!===============================================================================
subroutine read_ace_table(i_table, i_listing)
integer, intent(in) :: i_table ! index in nuclides/sab_tables
integer, intent(in) :: i_listing ! index in xs_listings
@ -258,9 +257,9 @@ contains
character(10) :: mat ! material identifier
character(70) :: comment ! comment for ACE table
character(MAX_FILE_LEN) :: filename ! path to ACE cross section library
type(Nuclide), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(XsListing), pointer :: listing => null()
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(XsListing), pointer :: listing
! determine path, record length, and location of table
listing => xs_listings(i_listing)
@ -406,8 +405,6 @@ contains
end select
deallocate(XSS)
if(associated(nuc)) nullify(nuc)
if(associated(sab)) nullify(sab)
end subroutine read_ace_table
@ -417,10 +414,8 @@ contains
!===============================================================================
subroutine read_esz(nuc, data_0K)
type(Nuclide), pointer :: nuc
logical :: data_0K ! are we reading 0K data?
type(Nuclide), intent(inout) :: nuc
logical, intent(in) :: data_0K ! are we reading 0K data?
integer :: NE ! number of energy points for total and elastic cross sections
integer :: i ! index in 0K elastic xs array for this nuclide
@ -507,8 +502,7 @@ contains
!===============================================================================
subroutine read_nu_data(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! loop index
integer :: JXS2 ! location for fission nu data
@ -524,7 +518,7 @@ contains
integer :: LOCC ! location of energy distributions for given MT
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist => null()
type(DistEnergy), pointer :: edist
JXS2 = JXS(2)
JXS24 = JXS(24)
@ -707,8 +701,7 @@ contains
!===============================================================================
subroutine read_reactions(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! loop indices
integer :: i_fission ! index in nuc % index_fission
@ -722,7 +715,6 @@ contains
integer :: IE ! reaction's starting index on energy grid
integer :: NE ! number of energies
integer :: NR ! number of interpolation regions
type(Reaction), pointer :: rxn => null()
type(ListInt) :: MTs
LMT = JXS(3)
@ -740,14 +732,15 @@ contains
! Store elastic scattering cross-section on reaction one -- note that the
! sigma array is not allocated or stored for elastic scattering since it is
! already stored in nuc % elastic
rxn => nuc % reactions(1)
rxn % MT = 2
rxn % Q_value = ZERO
rxn % multiplicity = 1
rxn % threshold = 1
rxn % scatter_in_cm = .true.
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
associate (rxn => nuc % reactions(1))
rxn%MT = 2
rxn%Q_value = ZERO
rxn%multiplicity = 1
rxn%threshold = 1
rxn%scatter_in_cm = .true.
rxn%has_angle_dist = .false.
rxn%has_energy_dist = .false.
end associate
! Add contribution of elastic scattering to total cross section
nuc % total = nuc % total + nuc % elastic
@ -760,123 +753,125 @@ contains
i_fission = 0
do i = 1, NMT
rxn => nuc % reactions(i+1)
associate (rxn => nuc % reactions(i+1))
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
! 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)
rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1)))
rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1)))
rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
! Read energy-dependent multiplicities
if (rxn % multiplicity > 100) then
! Set flag and allocate space for Tab1 to store yield
rxn % multiplicity_with_E = .true.
allocate(rxn % multiplicity_E)
! Read energy-dependent multiplicities
if (rxn % multiplicity > 100) then
! Set flag and allocate space for Tab1 to store yield
rxn % multiplicity_with_E = .true.
allocate(rxn % multiplicity_E)
XSS_index = JXS(11) + rxn % multiplicity - 101
NR = nint(XSS(XSS_index))
rxn % multiplicity_E % n_regions = NR
XSS_index = JXS(11) + rxn % multiplicity - 101
NR = nint(XSS(XSS_index))
rxn % multiplicity_E % n_regions = NR
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % multiplicity_E % nbt(NR))
allocate(rxn % multiplicity_E % int(NR))
end if
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % multiplicity_E % nbt(NR))
allocate(rxn % multiplicity_E % int(NR))
! read ENDF interpolation parameters
XSS_index = XSS_index + 1
if (NR > 0) then
rxn % multiplicity_E % nbt = get_int(NR)
rxn % multiplicity_E % int = get_int(NR)
end if
! allocate space for yield data
XSS_index = XSS_index + 2*NR
NE = nint(XSS(XSS_index))
rxn % multiplicity_E % n_pairs = NE
allocate(rxn % multiplicity_E % x(NE))
allocate(rxn % multiplicity_E % y(NE))
! read yield data
XSS_index = XSS_index + 1
rxn % multiplicity_E % x = get_real(NE)
rxn % multiplicity_E % y = get_real(NE)
end if
! read ENDF interpolation parameters
XSS_index = XSS_index + 1
if (NR > 0) then
rxn % multiplicity_E % nbt = get_int(NR)
rxn % multiplicity_E % int = get_int(NR)
end if
! read starting energy index
LOCA = int(XSS(LXS + i - 1))
IE = int(XSS(JXS7 + LOCA - 1))
rxn % threshold = IE
! allocate space for yield data
XSS_index = XSS_index + 2*NR
NE = nint(XSS(XSS_index))
rxn % multiplicity_E % n_pairs = NE
allocate(rxn % multiplicity_E % x(NE))
allocate(rxn % multiplicity_E % y(NE))
! read yield data
XSS_index = XSS_index + 1
rxn % multiplicity_E % x = get_real(NE)
rxn % multiplicity_E % y = get_real(NE)
end if
! read starting energy index
LOCA = int(XSS(LXS + i - 1))
IE = int(XSS(JXS7 + LOCA - 1))
rxn % threshold = IE
! read number of energies cross section values
NE = int(XSS(JXS7 + LOCA))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
! read number of energies cross section values
NE = int(XSS(JXS7 + LOCA))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
end associate
end do
! Create set of MT values
do i = 1, size(nuc % reactions)
call MTs % append(nuc % reactions(i) % MT)
call nuc%reaction_index%add_key(nuc%reactions(i)%MT, i)
end do
! Create total, absorption, and fission cross sections
do i = 2, size(nuc % reactions)
rxn => nuc % reactions(i)
IE = rxn % threshold
NE = size(rxn % sigma)
associate (rxn => nuc % reactions(i))
IE = rxn % threshold
NE = size(rxn % sigma)
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rxn % MT == N_LEVEL) cycle
if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rxn % MT == N_LEVEL) cycle
if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle
! Skip level cross sections if total is available
if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle
if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle
if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle
if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle
if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle
if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle
! Skip level cross sections if total is available
if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle
if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle
if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle
if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle
if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle
if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle
! Add contribution to total cross section
nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma
! Add contribution to total cross section
nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma
! Add contribution to absorption cross section
if (is_disappearance(rxn % MT)) then
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
end if
! Add contribution to absorption cross section
if (is_disappearance(rxn % MT)) then
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
end if
! Information about fission reactions
if (rxn % MT == N_FISSION) then
allocate(nuc % index_fission(1))
elseif (rxn % MT == N_F) then
allocate(nuc % index_fission(PARTIAL_FISSION_MAX))
nuc % has_partial_fission = .true.
end if
! Information about fission reactions
if (rxn % MT == N_FISSION) then
allocate(nuc % index_fission(1))
elseif (rxn % MT == N_F) then
allocate(nuc % index_fission(PARTIAL_FISSION_MAX))
nuc % has_partial_fission = .true.
end if
! Add contribution to fission cross section
if (is_fission(rxn % MT)) then
nuc % fissionable = .true.
nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
! Add contribution to fission cross section
if (is_fission(rxn % MT)) then
nuc % fissionable = .true.
nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
! Also need to add fission cross sections to absorption
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
! Also need to add fission cross sections to absorption
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to nuc % fission
if (rxn % MT == N_FISSION) deallocate(rxn % sigma)
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to nuc % fission
if (rxn % MT == N_FISSION) deallocate(rxn % sigma)
! Keep track of this reaction for easy searching later
i_fission = i_fission + 1
nuc % index_fission(i_fission) = i
nuc % n_fission = nuc % n_fission + 1
end if
! Keep track of this reaction for easy searching later
i_fission = i_fission + 1
nuc % index_fission(i_fission) = i
nuc % n_fission = nuc % n_fission + 1
end if
end associate
end do
! Clear MTs set
@ -890,8 +885,7 @@ contains
!===============================================================================
subroutine read_angular_dist(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
@ -902,7 +896,6 @@ contains
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
integer :: length ! length of data array to allocate
type(Reaction), pointer :: rxn => null()
JXS8 = JXS(8)
JXS9 = JXS(9)
@ -910,71 +903,72 @@ contains
! loop over all reactions with secondary neutrons -- NXS(5) does not include
! elastic scattering
do i = 1, NXS(5) + 1
rxn => nuc%reactions(i)
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 asssumed (in CM if TY < 0 and in LAB if TY > 0)
cycle
end if
rxn % has_angle_dist = .true.
! 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))
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! 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
! 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
end do
rxn % has_angle_dist = .true.
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! 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))
! 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)
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! 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
end do
! 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
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! 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)
! 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
end do
end associate
end do
end subroutine read_angular_dist
@ -985,29 +979,28 @@ contains
!===============================================================================
subroutine read_energy_dist(nuc)
type(Nuclide), pointer :: 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
type(Reaction), pointer :: rxn => null()
LED = JXS(10)
! Loop over all reactions
do i = 1, NXS(5)
rxn => nuc % reactions(i+1) ! skip over elastic scattering
rxn % has_energy_dist = .true.
associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering
rxn % has_energy_dist = .true.
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! allocate energy distribution
allocate(rxn % edist)
! allocate energy distribution
allocate(rxn % edist)
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
end associate
end do
end subroutine read_energy_dist
@ -1019,10 +1012,9 @@ contains
!===============================================================================
recursive subroutine get_energy_dist(edist, loc_law, delayed_n)
type(DistEnergy), pointer :: edist ! energy distribution
integer, intent(in) :: loc_law ! locator for data
logical, optional :: delayed_n ! is this for delayed neutrons?
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?
integer :: LDIS ! location of all energy distributions
integer :: LNW ! location of next energy distribution if multiple
@ -1102,7 +1094,6 @@ contains
!===============================================================================
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
@ -1146,7 +1137,7 @@ contains
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))
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
@ -1204,7 +1195,7 @@ contains
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))
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
@ -1234,7 +1225,7 @@ contains
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))
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
@ -1285,7 +1276,7 @@ contains
! 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))
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
@ -1301,8 +1292,7 @@ contains
!===============================================================================
subroutine read_unr_res(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: JXS23 ! location of URR data
integer :: lc ! locator
@ -1390,8 +1380,7 @@ contains
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide), pointer :: nuc
type(Nuclide), intent(inout) :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy
@ -1417,8 +1406,7 @@ contains
!===============================================================================
subroutine read_thermal_data(table)
type(SAlphaBeta), pointer :: table
type(SAlphaBeta), intent(inout) :: table
integer :: i ! index for incoming energies
integer :: j ! index for outgoing energies
@ -1600,7 +1588,7 @@ contains
do i = 1, n_nuclides_total
do j = 1, n_nuclides_total
if (nuclides(i) % zaid == nuclides(j) % zaid) then
call nuclides(i) % nuc_list % append(j)
call nuclides(i) % nuc_list % push_back(j)
end if
end do
end do

96
src/ace_header.F90
View file

@ -1,8 +1,9 @@
module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use endf_header, only: Tab1
use list_header, only: ListInt
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use stl_vector, only: VectorInt
implicit none
@ -17,10 +18,6 @@ module ace_header
integer, allocatable :: type(:) ! type of distribution
integer, allocatable :: location(:) ! location of each table
real(8), allocatable :: data(:) ! angular distribution data
! Type-Bound procedures
contains
procedure :: clear => distangle_clear ! Deallocates DistAngle
end type DistAngle
!===============================================================================
@ -51,7 +48,7 @@ module ace_header
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
@ -79,10 +76,6 @@ module ace_header
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
! Type-Bound procedures
contains
procedure :: clear => urrdata_clear ! Deallocates UrrData
end type UrrData
!===============================================================================
@ -99,7 +92,7 @@ module ace_header
real(8) :: kT ! temperature in MeV (k*T)
! Linked list of indices in nuclides array of instances of this same nuclide
type(ListInt) :: nuc_list
type(VectorInt) :: nuc_list
! Energy grid information
integer :: n_grid ! # of nuclide grid points
@ -153,7 +146,9 @@ module ace_header
! Reactions
integer :: n_reaction ! # of reactions
type(Reaction), pointer :: reactions(:) => null()
type(Reaction), allocatable :: reactions(:)
type(DictIntInt) :: reaction_index ! map MT values to index in reactions
! array; used at tally-time
! Type-Bound procedures
contains
@ -166,14 +161,12 @@ module ace_header
!===============================================================================
type Nuclide0K
character(10) :: nuclide ! name of nuclide, e.g. U-238
character(16) :: scheme = 'ares' ! target velocity sampling scheme
character(10) :: name ! name of nuclide, e.g. 92235.03c
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
real(8) :: E_min = 0.01e-6_8 ! lower cutoff energy for res scattering
real(8) :: E_max = 1000.0e-6_8 ! upper cutoff energy for res scattering
end type Nuclide0K
!===============================================================================
@ -265,7 +258,6 @@ module ace_header
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
real(8) :: kappa_fission ! microscopic energy-released from fission
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
@ -288,24 +280,10 @@ module ace_header
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
real(8) :: kappa_fission ! macroscopic energy-released from fission
end type MaterialMacroXS
contains
!===============================================================================
! DISTANGLE_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine distangle_clear(this)
class(DistAngle), intent(inout) :: this ! The DistAngle object to clear
if (allocated(this % energy)) &
deallocate(this % energy, this % type, this % location, this % data)
end subroutine distangle_clear
!===============================================================================
! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
!===============================================================================
@ -314,12 +292,6 @@ module ace_header
class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
! Clear p_valid
call this % p_valid % clear()
if (allocated(this % data)) &
deallocate(this % data)
if (associated(this % next)) then
! recursively clear this item
call this % next % clear()
@ -336,32 +308,13 @@ module ace_header
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (allocated(this % sigma)) deallocate(this % sigma)
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)
end if
call this % adist % clear()
end subroutine reaction_clear
!===============================================================================
! URRDATA_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine urrdata_clear(this)
class(UrrData), intent(inout) :: this ! The UrrData object to clear
if (allocated(this % energy)) &
deallocate(this % energy, this % prob)
end subroutine urrdata_clear
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide.
!===============================================================================
@ -372,31 +325,6 @@ module ace_header
integer :: i ! Loop counter
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &
& this % fission, this % nu_fission, this % absorption)
if (allocated(this % energy_0K)) &
deallocate(this % energy_0K)
if (allocated(this % elastic_0K)) &
deallocate(this % elastic_0K)
if (allocated(this % xs_cdf)) &
deallocate(this % xs_cdf)
if (allocated(this % heating)) &
deallocate(this % heating)
if (allocated(this % index_fission)) deallocate(this % index_fission)
if (allocated(this % nu_t_data)) deallocate(this % nu_t_data)
if (allocated(this % nu_p_data)) deallocate(this % nu_p_data)
if (allocated(this % nu_d_data)) deallocate(this % nu_d_data)
if (allocated(this % nu_d_precursor_data)) &
deallocate(this % nu_d_precursor_data)
if (associated(this % nu_d_edist)) then
do i = 1, size(this % nu_d_edist)
call this % nu_d_edist(i) % clear()
@ -405,18 +333,16 @@ module ace_header
end if
if (associated(this % urr_data)) then
call this % urr_data % clear()
deallocate(this % urr_data)
end if
if (associated(this % reactions)) then
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()
end do
deallocate(this % reactions)
end if
call this % nuc_list % clear()
call this % reaction_index % clear()
end subroutine nuclide_clear

11
src/constants.F90
View file

@ -8,7 +8,7 @@ module constants
! OpenMC major, minor, and release numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 7
integer, parameter :: VERSION_RELEASE = 0
integer, parameter :: VERSION_RELEASE = 1
! Revision numbers for binary files
integer, parameter :: REVISION_STATEPOINT = 14
@ -45,10 +45,11 @@ module constants
! Maximum number of words in a single line, length of line, and length of
! single word
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
integer, parameter :: REGION_SPEC_LEN = 1000
! Maximum number of external source spatial resamples to encounter before an
! error is thrown.

95
src/cross_section.F90
View file

@ -13,10 +13,6 @@ module cross_section
use search, only: binary_search
implicit none
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
@ -33,7 +29,8 @@ contains
integer :: i_nuclide ! index into nuclides array
integer :: i_sab ! index into sab_tables array
integer :: j ! index in mat % i_sab_nuclides
integer :: u ! index into logarithmic mapping array
integer :: i_grid ! index into logarithmic mapping array or material
! union grid
real(8) :: atom_density ! atom density of a nuclide
logical :: check_sab ! should we check for S(a,b) table?
type(Material), pointer :: mat ! current material
@ -44,7 +41,6 @@ contains
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % kappa_fission = ZERO
! Exit subroutine if material is void
if (p % material == MATERIAL_VOID) return
@ -52,11 +48,10 @@ contains
mat => materials(p % material)
! Find energy index on energy grid
u = 0
if (grid_method == GRID_MAT_UNION) then
call find_energy_index(p % E, p % material)
i_grid = find_energy_index(mat, p % E)
else if (grid_method == GRID_LOGARITHM) then
u = int(log(p % E/energy_min_neutron)/log_spacing)
i_grid = int(log(p % E/energy_min_neutron)/log_spacing)
end if
! Determine if this material has S(a,b) tables
@ -99,9 +94,9 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
end if
! ========================================================================
@ -129,10 +124,6 @@ contains
! Add contributions to material macroscopic nu-fission cross section
material_xs % nu_fission = material_xs % nu_fission + &
atom_density * micro_xs(i_nuclide) % nu_fission
! Add contributions to material macroscopic energy release from fission
material_xs % kappa_fission = material_xs % kappa_fission + &
atom_density * micro_xs(i_nuclide) % kappa_fission
end do
end subroutine calculate_xs
@ -142,18 +133,19 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, u)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_mat ! index into materials array
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer, intent(in) :: u ! index into logarithmic mapping array
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
real(8), intent(in) :: E ! energy
real(8) :: f ! interp factor on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
@ -165,7 +157,7 @@ contains
select case (grid_method)
case (GRID_MAT_UNION)
i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index)
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
@ -178,8 +170,8 @@ contains
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
i_low = nuc % grid_index(u)
i_high = nuc % grid_index(u + 1) + 1
i_low = nuc % grid_index(i_log_union)
i_high = nuc % grid_index(i_log_union + 1) + 1
! Perform binary search over reduced range
i_grid = binary_search(nuc % energy(i_low:i_high), &
@ -219,7 +211,6 @@ contains
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % kappa_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
@ -241,13 +232,6 @@ contains
! Calculate microscopic nuclide nu-fission cross section
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
i_grid) + f * nuc % nu_fission(i_grid+1)
! Calculate microscopic nuclide kappa-fission cross section
! The ENDF standard (ENDF-102) states that MT 18 stores
! the fission energy as the Q_value (fission(1))
micro_xs(i_nuclide) % kappa_fission = &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuclide) % fission
end if
! If there is S(a,b) data for this nuclide, we need to do a few
@ -382,7 +366,6 @@ contains
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer :: urr
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
micro_xs(i_nuclide) % use_ptable = .true.
@ -408,7 +391,7 @@ contains
! preserve correlation of temperature in probability tables
same_nuc = .false.
do i = 1, nuc % nuc_list % size()
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % data(i)) % last_E) then
same_nuc = .true.
same_nuc_idx = i
exit
@ -416,7 +399,7 @@ contains
end do
if (same_nuc) then
r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn
r = micro_xs(nuc % nuc_list % data(same_nuc_idx)) % last_prn
else
r = prn()
micro_xs(i_nuclide) % last_prn = r
@ -478,18 +461,17 @@ contains
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get pointer to inelastic scattering reaction
rxn => nuc % reactions(nuc % urr_inelastic)
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine inelastic scattering cross section
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
end associate
end if
! Multiply by smooth cross-section if needed
@ -526,38 +508,35 @@ contains
! energy
!===============================================================================
subroutine find_energy_index(E, i_mat)
real(8), intent(in) :: E ! energy of particle
integer, intent(in) :: i_mat ! material index
type(Material), pointer :: mat ! pointer to current material
mat => materials(i_mat)
pure function find_energy_index(mat, E) result(i)
type(Material), intent(in) :: mat ! pointer to current material
real(8), intent(in) :: E ! energy of particle
integer :: i ! energy grid index
! if the energy is outside of energy grid range, set to first or last
! index. Otherwise, do a binary search through the union energy grid.
if (E <= mat % e_grid(1)) then
union_grid_index = 1
i = 1
elseif (E > mat % e_grid(mat % n_grid)) then
union_grid_index = mat % n_grid - 1
i = mat % n_grid - 1
else
union_grid_index = binary_search(mat % e_grid, mat % n_grid, E)
i = binary_search(mat % e_grid, mat % n_grid, E)
end if
end subroutine find_energy_index
end function find_energy_index
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering
!===============================================================================
function elastic_xs_0K(E, nuc) result(xs_out)
pure function elastic_xs_0K(E, nuc) result(xs_out)
real(8), intent(in) :: E ! trial energy
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
type(Nuclide), pointer :: nuc ! target nuclide at temperature
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
real(8), intent(inout) :: E ! trial energy
real(8) :: xs_out ! 0K xs at trial energy
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
! Determine index on nuclide energy grid
if (E < nuc % energy_0K(1)) then

2
src/endf.F90
View file

@ -11,7 +11,7 @@ contains
! REACTION_NAME gives the name of the reaction for a given MT value
!===============================================================================
function reaction_name(MT) result(string)
pure function reaction_name(MT) result(string)
integer, intent(in) :: MT
character(20) :: string

22
src/endf_header.F90
View file

@ -13,28 +13,6 @@ 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
! Type-Bound procedures
contains
procedure :: clear => tab1_clear ! deallocates a Tab1 Object.
end type Tab1
contains
!===============================================================================
! TAB1_CLEAR deallocates the items in Tab1
!===============================================================================
subroutine tab1_clear(this)
class(Tab1), intent(inout) :: this ! The Tab1 to clear
if (allocated(this % nbt)) &
deallocate(this % nbt, this % int)
if (allocated(this % x)) &
deallocate(this % x, this % y)
end subroutine tab1_clear
end module endf_header

28
src/fission.F90
View file

@ -15,18 +15,17 @@ contains
! given nuclide and incoming neutron energy
!===============================================================================
function nu_total(nuc, E) result(nu)
type(Nuclide), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
pure function nu_total(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_t_type == NU_NONE) then
call fatal_error("No neutron emission data for table: " // nuc % name)
nu = ERROR_REAL
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
@ -49,11 +48,10 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_prompt(nuc, E) result(nu)
type(Nuclide), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
pure function nu_prompt(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
@ -87,10 +85,9 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_delayed(nuc, E) result(nu)
pure function nu_delayed(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
if (nuc % nu_d_type == NU_NONE) then
@ -111,8 +108,7 @@ contains
! a given nuclide and incoming neutron energy in a given delayed group.
!===============================================================================
function yield_delayed(nuc, E, g) result(yield)
pure function yield_delayed(nuc, E, g) result(yield)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: yield ! delayed neutron precursor yield

126
src/geometry.F90
View file

@ -9,6 +9,7 @@ module geometry
use particle_header, only: LocalCoord, Particle
use particle_restart_write, only: write_particle_restart
use surface_header
use stl_vector, only: VectorInt
use string, only: to_str
use tally, only: score_surface_current
@ -119,6 +120,7 @@ contains
stack(i_stack) = (actual_sense .eqv. (token > 0))
end if
end select
end do
if (i_stack == 1) then
@ -597,8 +599,9 @@ contains
real(8) :: d_lat ! distance to lattice boundary
real(8) :: d_surf ! distance to surface
real(8) :: x0,y0,z0 ! coefficients for surface
real(8) :: xyz_cross(3) ! coordinates at projected surface crossing
logical :: coincident ! is particle on surface?
type(Cell), pointer :: cl
type(Cell), pointer :: c
class(Surface), pointer :: surf
class(Lattice), pointer :: lat
@ -614,7 +617,7 @@ contains
LEVEL_LOOP: do j = 1, p % n_coord
! get pointer to cell on this level
cl => cells(p % coord(j) % cell)
c => cells(p % coord(j) % cell)
! copy directional cosines
u = p % coord(j) % uvw(1)
@ -624,8 +627,8 @@ contains
! =======================================================================
! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
SURFACE_LOOP: do i = 1, size(cl%region)
index_surf = cl%region(i)
SURFACE_LOOP: do i = 1, size(c % region)
index_surf = c % region(i)
coincident = (index_surf == p % surface)
! ignore this token if it corresponds to an operator rather than a
@ -634,14 +637,14 @@ contains
if (index_surf >= OP_UNION) cycle
! Calculate distance to surface
surf => surfaces(index_surf)%obj
d = surf%distance(p%coord(j)%xyz, p%coord(j)%uvw, coincident)
surf => surfaces(index_surf) % obj
d = surf % distance(p % coord(j) % xyz, p % coord(j) % uvw, coincident)
! Check if calculated distance is new minimum
if (d < d_surf) then
if (abs(d - d_surf)/d_surf >= FP_PRECISION) then
d_surf = d
level_surf_cross = -cl % region(i)
level_surf_cross = -c % region(i)
end if
end if
end do SURFACE_LOOP
@ -847,14 +850,31 @@ contains
if (d_surf < d_lat) then
if ((dist - d_surf)/dist >= FP_REL_PRECISION) then
dist = d_surf
surface_crossed = level_surf_cross
! If the cell is not simple, it is possible that both the negative and
! positive half-space were given in the region specification. Thus, we
! have to explicitly check which half-space the particle would be
! traveling into if the surface is crossed
if (.not. c % simple) then
xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw
surf => surfaces(abs(level_surf_cross)) % obj
if (dot_product(p % coord(j) % uvw, &
surf % normal(xyz_cross)) > ZERO) then
surface_crossed = abs(level_surf_cross)
else
surface_crossed = -abs(level_surf_cross)
end if
else
surface_crossed = level_surf_cross
end if
lattice_translation(:) = [0, 0, 0]
next_level = j
end if
else
if ((dist - d_lat)/dist >= FP_REL_PRECISION) then
dist = d_lat
surface_crossed = None
surface_crossed = NONE
lattice_translation(:) = level_lat_trans
next_level = j
end if
@ -871,81 +891,51 @@ contains
subroutine neighbor_lists()
integer :: i ! index in cells/surfaces array
integer :: j ! index of surface in cell
integer :: i_surface ! index in count arrays
integer, allocatable :: count_positive(:) ! # of cells on positive side
integer, allocatable :: count_negative(:) ! # of cells on negative side
logical :: positive ! positive side specified in surface list
type(Cell), pointer :: c
integer :: i ! index in cells/surfaces array
integer :: j ! index in region specification
integer :: k ! surface half-space spec
integer :: n ! size of vector
type(VectorInt), allocatable :: neighbor_pos(:)
type(VectorInt), allocatable :: neighbor_neg(:)
call write_message("Building neighboring cells lists for each surface...", &
&4)
4)
allocate(count_positive(n_surfaces))
allocate(count_negative(n_surfaces))
count_positive = 0
count_negative = 0
allocate(neighbor_pos(n_surfaces))
allocate(neighbor_neg(n_surfaces))
do i = 1, n_cells
c => cells(i)
do j = 1, size(cells(i)%region)
! Get token from region specification and skip any tokens that
! correspond to operators rather than regions
k = cells(i)%region(j)
if (abs(k) >= OP_UNION) cycle
! loop over each region specification
do j = 1, size(c%region)
i_surface = c % region(j)
positive = (i_surface > 0)
! Skip any tokens that correspond to operators rather than regions
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1
! Add this cell ID to neighbor list for k-th surface
if (k > 0) then
call neighbor_pos(abs(k))%push_back(i)
else
count_negative(i_surface) = count_negative(i_surface) + 1
call neighbor_neg(abs(k))%push_back(i)
end if
end do
end do
! allocate neighbor lists for each surface
do i = 1, n_surfaces
if (count_positive(i) > 0) then
allocate(surfaces(i)%obj%neighbor_pos(count_positive(i)))
! Copy positive neighbors to Surface instance
n = neighbor_pos(i)%size()
if (n > 0) then
allocate(surfaces(i)%obj%neighbor_pos(n))
surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n)
end if
if (count_negative(i) > 0) then
allocate(surfaces(i)%obj%neighbor_neg(count_negative(i)))
! Copy negative neighbors to Surface instance
n = neighbor_neg(i)%size()
if (n > 0) then
allocate(surfaces(i)%obj%neighbor_neg(n))
surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n)
end if
end do
count_positive = 0
count_negative = 0
! loop over all cells
do i = 1, n_cells
c => cells(i)
! loop through the region specification
do j = 1, size(c%region)
i_surface = c % region(j)
positive = (i_surface > 0)
! Skip any tokens that correspond to operators rather than regions
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1
surfaces(i_surface)%obj%neighbor_pos(count_positive(i_surface)) = i
else
count_negative(i_surface) = count_negative(i_surface) + 1
surfaces(i_surface)%obj%neighbor_neg(count_negative(i_surface)) = i
end if
end do
end do
deallocate(count_positive)
deallocate(count_negative)
end subroutine neighbor_lists
!===============================================================================

16
src/global.F90
View file

@ -436,7 +436,12 @@ contains
do i = 1, size(nuclides)
call nuclides(i) % clear()
end do
deallocate(nuclides)
! WARNING: The following statement should work but doesn't under gfortran
! 4.6 because of a bug. Technically, commenting this out leaves a memory
! leak.
! deallocate(nuclides)
end if
if (allocated(nuclides_0K)) then
@ -463,14 +468,7 @@ contains
! Deallocate tally-related arrays
if (allocated(global_tallies)) deallocate(global_tallies)
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) then
! First call the clear routines
do i = 1, size(tallies)
call tallies(i) % clear()
end do
! Now deallocate the tally array
deallocate(tallies)
end if
if (allocated(tallies)) deallocate(tallies)
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(tally_maps)) deallocate(tally_maps)

14
src/hdf5_interface.F90
View file

@ -1483,7 +1483,7 @@ contains
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: filetype
integer(HID_T) :: memtype
integer(HSIZE_T) :: n
integer(SIZE_T) :: n
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1544,8 +1544,8 @@ contains
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: filetype
integer(HID_T) :: memtype
integer(HSIZE_T) :: size
integer(HSIZE_T) :: n
integer(SIZE_T) :: size
integer(SIZE_T) :: n
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1628,7 +1628,7 @@ contains
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: filetype
integer(HID_T) :: memtype
integer(HSIZE_T) :: n
integer(SIZE_T) :: n
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O
@ -1644,7 +1644,7 @@ contains
! Create datatype in memory based on Fortran character
call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err)
call h5tset_size_f(memtype, int(len(buffer(1)), HSIZE_T), hdf5_err)
call h5tset_size_f(memtype, int(len(buffer(1)), SIZE_T), hdf5_err)
! Create dataspace/dataset
call h5screate_simple_f(1, dims, dspace, hdf5_err)
@ -1706,8 +1706,8 @@ contains
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: filetype
integer(HID_T) :: memtype
integer(HSIZE_T) :: size
integer(HSIZE_T) :: n
integer(SIZE_T) :: size
integer(SIZE_T) :: n
type(c_ptr) :: f_ptr
! Set up collective vs. independent I/O

2
src/input_xml.F90
View file

@ -994,7 +994,7 @@ contains
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
character(1000) :: region_spec
character(REGION_SPEC_LEN) :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
class(Lattice), pointer :: lat

32
src/interpolation.F90
View file

@ -21,7 +21,7 @@ contains
! tabulated x's and y's.
!===============================================================================
function interpolate_tab1_array(data, x, loc_start) result(y)
pure function interpolate_tab1_array(data, x, loc_start) result(y)
real(8), intent(in) :: data(:) ! array of data
real(8), intent(in) :: x ! x value to find y at
@ -106,18 +106,16 @@ contains
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = (1 - r)*y0 + r*y1
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = (1 - r)*y0 + r*y1
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = exp((1-r)*log(y0) + r*log(y1))
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
call fatal_error("Unsupported interpolation scheme: " // to_str(interp))
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_array
@ -129,7 +127,7 @@ contains
! tabulated x's and y's.
!===============================================================================
function interpolate_tab1_object(obj, x) result(y)
pure function interpolate_tab1_object(obj, x) result(y)
type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function
real(8), intent(in) :: x ! x value to find y at
@ -191,18 +189,16 @@ contains
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = (1 - r)*y0 + r*y1
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = (1 - r)*y0 + r*y1
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = exp((1-r)*log(y0) + r*log(y1))
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
call fatal_error("Unsupported interpolation scheme: " // to_str(interp))
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_object

6
src/math.F90
View file

@ -12,7 +12,7 @@ contains
! distribution with a specified probability level
!===============================================================================
function normal_percentile(p) result(z)
elemental function normal_percentile(p) result(z)
real(8), intent(in) :: p ! probability level
real(8) :: z ! corresponding z-value
@ -71,7 +71,7 @@ contains
! specified probability level and number of degrees of freedom
!===============================================================================
function t_percentile(p, df) result(t)
elemental function t_percentile(p, df) result(t)
real(8), intent(in) :: p ! probability level
integer, intent(in) :: df ! degrees of freedom
@ -123,7 +123,7 @@ contains
! the return value will be 1.0.
!===============================================================================
pure function calc_pn(n,x) result(pnx)
elemental function calc_pn(n,x) result(pnx)
integer, intent(in) :: n ! Legendre order requested
real(8), intent(in) :: x ! Independent variable the Legendre is to be

38
src/mesh.F90
View file

@ -18,9 +18,8 @@ contains
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
subroutine get_mesh_bin(m, xyz, bin)
type(RegularMesh), pointer :: m ! mesh pointer
pure subroutine get_mesh_bin(m, xyz, bin)
type(RegularMesh), intent(in) :: m ! mesh pointer
real(8), intent(in) :: xyz(:) ! coordinates
integer, intent(out) :: bin ! tally bin
@ -71,9 +70,8 @@ contains
! GET_MESH_INDICES determines the indices of a particle in a structured mesh
!===============================================================================
subroutine get_mesh_indices(m, xyz, ijk, in_mesh)
type(RegularMesh), pointer :: m
pure subroutine get_mesh_indices(m, xyz, ijk, in_mesh)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz(:) ! coordinates to check
integer, intent(out) :: ijk(:) ! indices in mesh
logical, intent(out) :: in_mesh ! were given coords in mesh?
@ -96,11 +94,10 @@ contains
! use in a TallyObject results array
!===============================================================================
function mesh_indices_to_bin(m, ijk, surface_current) result(bin)
type(RegularMesh), pointer :: m
pure function mesh_indices_to_bin(m, ijk, surface_current) result(bin)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: ijk(:)
logical, optional :: surface_current
logical, intent(in), optional :: surface_current
integer :: bin
integer :: n_y ! number of mesh cells in y direction
@ -130,9 +127,8 @@ contains
! (i,j) or (i,j,k) indices
!===============================================================================
subroutine bin_to_mesh_indices(m, bin, ijk)
type(RegularMesh), pointer :: m
pure subroutine bin_to_mesh_indices(m, bin, ijk)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: bin
integer, intent(out) :: ijk(:)
@ -167,9 +163,9 @@ contains
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
! cell and energy group
real(8), optional :: energies(:) ! energy grid to search
integer(8), optional :: size_bank ! # of bank sites (on each proc)
logical, optional :: sites_outside ! were there sites outside mesh?
real(8), intent(in), optional :: energies(:) ! energy grid to search
integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc)
logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?
integer :: i ! loop index for local fission sites
integer :: n_sites ! size of bank array
@ -262,9 +258,8 @@ contains
! track will score to a mesh tally.
!===============================================================================
function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), pointer :: m
pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(2)
real(8), intent(in) :: xyz1(2)
logical :: intersects
@ -328,9 +323,8 @@ contains
end function mesh_intersects_2d
function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), pointer :: m
pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(3)
real(8), intent(in) :: xyz1(3)
logical :: intersects

99
src/output.F90
View file

@ -100,10 +100,9 @@ contains
!===============================================================================
subroutine header(msg, unit, level)
character(*), intent(in) :: msg ! header message
integer, optional :: unit ! unit to write to
integer, optional :: level ! specified header level
integer, intent(in), optional :: unit ! unit to write to
integer, intent(in), optional :: level ! specified header level
integer :: n ! number of = signs on left
integer :: m ! number of = signs on right
@ -195,9 +194,8 @@ contains
!===============================================================================
subroutine write_message(message, level)
character(*) :: message
integer, optional :: level ! verbosity level
character(*), intent(in) :: message ! message to write
integer, intent(in), optional :: level ! verbosity level
integer :: i_start ! starting position
integer :: i_end ! ending position
@ -250,7 +248,6 @@ contains
!===============================================================================
subroutine print_particle(p)
type(Particle), intent(in) :: p
integer :: i ! index for coordinate levels
@ -320,9 +317,8 @@ contains
!===============================================================================
subroutine print_nuclide(nuc, unit)
type(Nuclide), pointer :: nuc
integer, optional :: unit
type(Nuclide), intent(in) :: nuc
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
@ -334,8 +330,7 @@ contains
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(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
type(UrrData), pointer :: urr
! set default unit for writing information
if (present(unit)) then
@ -363,32 +358,32 @@ contains
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
do i = 1, nuc % n_reaction
rxn => nuc % reactions(i)
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 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
! 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)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
! 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
! 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
! Add memory required for summary reactions (total, absorption, fission,
@ -438,9 +433,8 @@ contains
!===============================================================================
subroutine print_sab_table(sab, unit)
type(SAlphaBeta), pointer :: sab
integer, optional :: unit
type(SAlphaBeta), intent(in) :: sab
integer, intent(in), optional :: unit
integer :: size_sab ! memory used by S(a,b) table
integer :: unit_ ! unit to write to
@ -526,8 +520,8 @@ contains
integer :: i ! loop index
integer :: unit_xs ! cross_sections.out file unit
character(MAX_FILE_LEN) :: path ! path of summary file
type(Nuclide), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
! Create filename for log file
path = trim(path_output) // "cross_sections.out"
@ -681,7 +675,7 @@ contains
subroutine print_plot()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl => null()
type(ObjectPlot), pointer :: pl
! Display header for plotting
call header("PLOTTING SUMMARY")
@ -1201,7 +1195,7 @@ contains
!===============================================================================
subroutine write_surface_current(t, unit_tally)
type(TallyObject), pointer :: t
type(TallyObject), intent(in) :: t
integer, intent(in) :: unit_tally
integer :: i ! mesh index for x
@ -1373,10 +1367,9 @@ contains
!===============================================================================
function get_label(t, i_filter) result(label)
type(TallyObject), pointer :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(100) :: label ! user-specified identifier
type(TallyObject), intent(in) :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(100) :: label ! user-specified identifier
integer :: i ! index in cells/surfaces/etc array
integer :: bin
@ -1440,12 +1433,12 @@ contains
recursive subroutine find_offset(map, goal, univ, final, offset, path)
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
type(Universe), pointer, intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(100) :: path ! Path to offset
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
@ -1457,7 +1450,7 @@ contains
integer :: temp_offset ! Looped sum of offsets
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
type(Cell), pointer:: c ! Pointer to current cell
type(Cell), pointer :: c ! Pointer to current cell
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice

86
src/physics.F90
View file

@ -185,7 +185,6 @@ contains
!===============================================================================
subroutine sample_fission(i_nuclide, i_reaction)
integer, intent(in) :: i_nuclide ! index in nuclides array
integer, intent(out) :: i_reaction ! index in nuc % reactions array
@ -195,7 +194,6 @@ contains
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -220,14 +218,15 @@ contains
FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
i_reaction = nuc % index_fission(i)
rxn => nuc % reactions(i_reaction)
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
associate (rxn => nuc % reactions(i_reaction))
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
end associate
! Create fission bank sites if fission occurs
if (prob > cutoff) exit FISSION_REACTION_LOOP
@ -312,11 +311,10 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
type(Nuclide), pointer :: nuc
! copy incoming direction
uvw_old(:) = p % coord(1) % uvw
@ -343,11 +341,8 @@ contains
p % E, p % coord(1) % uvw, p % mu)
else
! get pointer to elastic scattering reaction
rxn => nuc % reactions(1)
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, rxn, &
call elastic_scatter(i_nuclide, nuc % reactions(1), &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
end if
@ -370,28 +365,28 @@ contains
&// trim(nuc % name))
end if
rxn => nuc % reactions(i)
associate (rxn => nuc % reactions(i))
! Skip fission reactions
if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF &
.or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle
! Skip fission reactions
if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF &
.or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle
! some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle
! some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
end associate
end do
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, rxn, p)
p % event_MT = rxn % MT
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc%reactions(i)%MT
end if
@ -420,9 +415,8 @@ contains
!===============================================================================
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt)
integer, intent(in) :: i_nuclide
type(Reaction), pointer :: rxn
type(Reaction), intent(in) :: rxn
real(8), intent(inout) :: E
real(8), intent(inout) :: uvw(3)
real(8), intent(out) :: mu_lab
@ -759,9 +753,7 @@ contains
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
type(Nuclide), pointer :: nuc ! target nuclide at temperature T
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T
real(8), intent(out) :: v_target(3) ! target velocity
real(8), intent(in) :: v_neut(3) ! neutron velocity
real(8), intent(in) :: E ! particle energy
@ -1006,8 +998,7 @@ contains
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
type(Nuclide), pointer :: nuc ! target nuclide at temperature
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature
real(8), intent(out) :: v_target(3)
real(8), intent(in) :: E
real(8), intent(in) :: uvw(3)
@ -1080,7 +1071,6 @@ contains
!===============================================================================
subroutine create_fission_sites(p, i_nuclide, i_reaction)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_reaction
@ -1095,11 +1085,9 @@ contains
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointers
nuc => nuclides(i_nuclide)
rxn => nuc % reactions(i_reaction)
! TODO: Heat generation from fission
@ -1170,7 +1158,8 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
fission_bank(i) % E = sample_fission_energy(nuc, rxn, p)
fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(&
i_reaction), p)
! Set the delayed group of the neutron
fission_bank(i) % delayed_group = p % delayed_group
@ -1197,8 +1186,8 @@ contains
function sample_fission_energy(nuc, rxn, p) result(E_out)
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing energy of fission neutron
@ -1323,8 +1312,8 @@ contains
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, p)
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
integer :: i ! loop index
@ -1409,8 +1398,7 @@ contains
!===============================================================================
function sample_angle(rxn, E) result(mu)
type(Reaction), pointer :: rxn ! reaction
type(Reaction), intent(in) :: rxn ! reaction
real(8), intent(in) :: E ! incoming energy
real(8) :: xi ! random number on [0,1)
@ -1536,7 +1524,6 @@ contains
!===============================================================================
function rotate_angle(uvw0, mu) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8) :: uvw(3) ! rotated directional cosine
@ -1585,8 +1572,7 @@ contains
!===============================================================================
recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
type(DistEnergy), pointer :: edist
type(DistEnergy), intent(in) :: edist
real(8), intent(in) :: E_in ! incoming energy of neutron
real(8), intent(out) :: E_out ! outgoing energy
real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle

27
src/search.F90
View file

@ -18,7 +18,7 @@ contains
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
function binary_search_real(array, n, val) result(array_index)
pure function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
@ -33,7 +33,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -49,8 +50,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do
@ -58,7 +59,7 @@ contains
end function binary_search_real
function binary_search_int4(array, n, val) result(array_index)
pure function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
@ -73,7 +74,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -89,8 +91,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do
@ -98,7 +100,7 @@ contains
end function binary_search_int4
function binary_search_int8(array, n, val) result(array_index)
pure function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
@ -113,7 +115,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -129,8 +132,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do

3
src/source.F90
View file

@ -96,8 +96,7 @@ contains
!===============================================================================
subroutine sample_external_source(site)
type(Bank), pointer :: site ! source site
type(Bank), intent(inout) :: site ! source site
integer :: i ! dummy loop index
real(8) :: r(3) ! sampled coordinates

2
src/state_point.F90
View file

@ -900,7 +900,6 @@ contains
integer(HSIZE_T) :: dims(1)
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer :: data_xfer_mode
integer(HID_T) :: plist ! property list
#else
integer :: i
@ -1019,7 +1018,6 @@ contains
integer(HSIZE_T) :: offset(1) ! offset of data
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer :: data_xfer_mode
integer(HID_T) :: plist ! property list
#endif

82
src/string.F90
View file

@ -25,7 +25,6 @@ contains
!===============================================================================
subroutine split_string(string, words, n)
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
@ -166,7 +165,7 @@ contains
! string = concatenated string
!===============================================================================
function concatenate(words, n_words) result(string)
pure function concatenate(words, n_words) result(string)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words)
@ -186,8 +185,7 @@ contains
! TO_LOWER converts a string to all lower case characters
!===============================================================================
function to_lower(word) result(word_lower)
pure function to_lower(word) result(word_lower)
character(*), intent(in) :: word
character(len=len(word)) :: word_lower
@ -209,8 +207,7 @@ contains
! TO_UPPER converts a string to all upper case characters
!===============================================================================
function to_upper(word) result(word_upper)
pure function to_upper(word) result(word_upper)
character(*), intent(in) :: word
character(len=len(word)) :: word_upper
@ -234,39 +231,38 @@ contains
! integers.
!===============================================================================
function zero_padded(num, n_digits) result(str)
integer, intent(in) :: num
integer, intent(in) :: n_digits
character(11) :: str
function zero_padded(num, n_digits) result(str)
integer, intent(in) :: num
integer, intent(in) :: n_digits
character(11) :: str
character(8) :: zp_form
character(8) :: zp_form
! Make sure n_digits is reasonable. 10 digits is the maximum needed for the
! largest integer(4).
if (n_digits > 10) then
call fatal_error('zero_padded called with an unreasonably large &
&n_digits (>10)')
end if
! Make sure n_digits is reasonable. 10 digits is the maximum needed for the
! largest integer(4).
if (n_digits > 10) then
call fatal_error('zero_padded called with an unreasonably large &
&n_digits (>10)')
end if
! Write a format string of the form '(In.m)' where n is the max width and
! m is the min width. If a sign is present, then n must be one greater
! than m.
if (num < 0) then
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits
else
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits
end if
! Write a format string of the form '(In.m)' where n is the max width and
! m is the min width. If a sign is present, then n must be one greater
! than m.
if (num < 0) then
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits+1, n_digits
else
write(zp_form, '("(I", I0, ".", I0, ")")') n_digits, n_digits
end if
! Format the number.
write(str, zp_form) num
end function zero_padded
! Format the number.
write(str, zp_form) num
end function zero_padded
!===============================================================================
! IS_NUMBER determines whether a string of characters is all 0-9 characters
!===============================================================================
function is_number(word) result(number)
pure function is_number(word) result(number)
character(*), intent(in) :: word
logical :: number
@ -286,10 +282,9 @@ end function zero_padded
! sequence of characters
!===============================================================================
logical function starts_with(str, seq)
character(*) :: str ! string to check
character(*) :: seq ! sequence of characters
pure logical function starts_with(str, seq)
character(*), intent(in) :: str ! string to check
character(*), intent(in) :: seq ! sequence of characters
integer :: i
integer :: i_start
@ -321,10 +316,9 @@ end function zero_padded
! of characters
!===============================================================================
logical function ends_with(str, seq)
character(*) :: str ! string to check
character(*) :: seq ! sequence of characters
pure logical function ends_with(str, seq)
character(*), intent(in) :: str ! string to check
character(*), intent(in) :: seq ! sequence of characters
integer :: i_start
integer :: str_len
@ -350,7 +344,7 @@ end function zero_padded
! integer.
!===============================================================================
function count_digits(num) result(n_digits)
pure function count_digits(num) result(n_digits)
integer, intent(in) :: num
integer :: n_digits
@ -368,7 +362,7 @@ end function zero_padded
! INT4_TO_STR converts an integer(4) to a string.
!===============================================================================
function int4_to_str(num) result(str)
pure function int4_to_str(num) result(str)
integer, intent(in) :: num
character(11) :: str
@ -382,7 +376,7 @@ end function zero_padded
! INT8_TO_STR converts an integer(8) to a string.
!===============================================================================
function int8_to_str(num) result(str)
pure function int8_to_str(num) result(str)
integer(8), intent(in) :: num
character(21) :: str
@ -396,7 +390,7 @@ end function zero_padded
! STR_TO_INT converts a string to an integer.
!===============================================================================
function str_to_int(str) result(num)
pure function str_to_int(str) result(num)
character(*), intent(in) :: str
integer(8) :: num
@ -421,7 +415,7 @@ end function zero_padded
! STR_TO_REAL converts an arbitrary string to a real(8)
!===============================================================================
function str_to_real(string) result(num)
pure function str_to_real(string) result(num)
character(*), intent(in) :: string
real(8) :: num
@ -440,7 +434,7 @@ end function zero_padded
! are used.
!===============================================================================
function real_to_str(num, sig_digits) result(string)
pure function real_to_str(num, sig_digits) result(string)
real(8), intent(in) :: num ! number to convert
integer, optional, intent(in) :: sig_digits ! # of significant digits

2
src/summary.F90
View file

@ -113,7 +113,7 @@ contains
integer(HID_T) :: universes_group, univ_group
integer(HID_T) :: lattices_group, lattice_group
real(8), allocatable :: coeffs(:)
character(MAX_LINE_LEN) :: region_spec
character(REGION_SPEC_LEN) :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
type(Universe), pointer :: u

280
src/tally.F90
View file

@ -37,13 +37,13 @@ contains
subroutine score_general(p, t, start_index, filter_index, i_nuclide, &
atom_density, flux)
type(Particle), intent(in) :: p
type(TallyObject), pointer, intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
integer :: i ! loop index for scoring bins
integer :: l ! loop index for nuclides in material
@ -67,9 +67,6 @@ contains
real(8) :: uvw(3) ! particle direction
real(8) :: E ! particle energy
logical :: scoring_diff_nuclide
type(Material), pointer :: mat
type(Reaction), pointer :: rxn
type(Nuclide), pointer :: nuc
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -221,24 +218,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -258,24 +251,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -295,24 +284,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -467,9 +452,6 @@ contains
! delayed-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
! Get the event nuclide
nuc => nuclides(p % event_nuclide)
! Check if the delayed group filter is present
if (dg_filter > 0) then
@ -481,11 +463,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(p % event_nuclide), E, d)
! Compute the score and tally to bin
score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuc, E) / &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
@ -495,7 +477,7 @@ contains
! by multiplying the absorbed weight by the fraction of the
! delayed-nu-fission xs to the absorption xs
score = p % absorb_wgt * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuc, E) / &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
end if
end if
@ -536,9 +518,6 @@ contains
! Check if tally is on a single nuclide
if (i_nuclide > 0) then
! Get the nuclide of interest
nuc => nuclides(i_nuclide)
! Check if the delayed group filter is present
if (dg_filter > 0) then
@ -549,11 +528,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(i_nuclide), E, d)
! Compute the score and tally to bin
score = micro_xs(i_nuclide) % fission * yield &
* nu_delayed(nuc, E) * atom_density * flux
* nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
@ -561,27 +540,24 @@ contains
! If the delayed group filter is not present, compute the score
! by multiplying the delayed-nu-fission macro xs by the flux
score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, E)&
* atom_density * flux
score = micro_xs(i_nuclide) % fission * &
nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
end if
! Tally is on total nuclides
else
! Get pointer to current material
mat => materials(p % material)
! Check if the delayed group filter is present
if (dg_filter > 0) then
! Loop over all nuclides in the current material
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
i_nuc = materials(p % material) % nuclide(l)
! Loop over all delayed group bins and tally to them individually
do d_bin = 1, t % filters(dg_filter) % n_bins
@ -589,15 +565,12 @@ contains
! Get the delayed group for this bin
d = t % filters(dg_filter) % int_bins(d_bin)
! Get the current nuclide
nuc => nuclides(i_nuc)
! Get the yield for the desired nuclide and delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(i_nuc), E, d)
! Compute the score and tally to bin
score = micro_xs(i_nuc) % fission * yield &
* nu_delayed(nuc, E) * atom_density_ * flux
* nu_delayed(nuclides(i_nuc), E) * atom_density_ * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
end do
@ -607,13 +580,13 @@ contains
score = ZERO
! Loop over all nuclides in the current material
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
i_nuc = materials(p % material) % nuclide(l)
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission &
@ -625,38 +598,67 @@ contains
case (SCORE_KAPPA_FISSION)
! Determine kappa-fission cross section on the fly. The ENDF standard
! (ENDF-102) states that MT 18 stores the fission energy as the Q_value
! (fission(1))
score = ZERO
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scale by kappa-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
score = p % absorb_wgt * &
micro_xs(p % event_nuclide) % kappa_fission / &
micro_xs(p % event_nuclide) % absorption
else
score = ZERO
end if
associate (nuc => nuclides(p%event_nuclide))
if (micro_xs(p%event_nuclide)%absorption > ZERO .and. &
nuc%fissionable) then
score = p%absorb_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
score = p % last_wgt * &
micro_xs(p % event_nuclide) % kappa_fission / &
micro_xs(p % event_nuclide) % absorption
associate (nuc => nuclides(p%event_nuclide))
if (nuc%fissionable) then
score = p%last_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
end if
end associate
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux
associate (nuc => nuclides(i_nuclide))
if (nuc%fissionable) then
score = nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuclide)%fission * atom_density * flux
end if
end associate
else
score = material_xs % kappa_fission * flux
do l = 1, materials(p%material)%n_nuclides
! Determine atom density and index of nuclide
atom_density_ = materials(p%material)%atom_density(l)
i_nuc = materials(p%material)%nuclide(l)
! If nuclide is fissionable, accumulate kappa fission
associate(nuc => nuclides(i_nuc))
if (nuc % fissionable) then
score = score + nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuc)%fission * atom_density_ * flux
end if
end associate
end do
end if
end if
case (SCORE_EVENTS)
! Simply count number of scoring events
score = ONE
@ -699,14 +701,10 @@ contains
score = ZERO
if (i_nuclide > 0) then
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide instance
! of the form {MT: i_reaction, ...}
REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuclide) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuclide)%reaction_index%get_key(score_bin)
associate (rxn => nuclides(i_nuclide) % reactions(m))
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuclide) % index_grid
@ -716,26 +714,20 @@ contains
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit REACTION_LOOP
end if
end do REACTION_LOOP
end associate
end if
else
! Get pointer to current material
mat => materials(p % material)
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide
! instance of the form {MT: i_reaction, ...}
do m = 1, nuclides(i_nuc) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuc) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
i_nuc = materials(p % material) % nuclide(l)
if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuc)%reaction_index%get_key(score_bin)
associate (rxn => nuclides(i_nuc) % reactions(m))
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuc) % index_grid
@ -745,9 +737,8 @@ contains
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density_ * flux
end if
exit
end if
end do
end associate
end if
end do
end if
@ -771,8 +762,8 @@ contains
case (ESTIMATOR_COLLISION)
if (t % deriv % dep_var == DIFF_NUCLIDE_DENSITY) then
mat => materials(p % material)
scoring_diff_nuclide = (mat % id == t % deriv % diff_material) &
scoring_diff_nuclide = &
(materials(p % material) % id == t % deriv % diff_material) &
.and. (i_nuclide == t % deriv % diff_nuclide)
select case (score_bin)
case (SCORE_TOTAL)
@ -787,10 +778,6 @@ contains
case (SCORE_NU_FISSION)
scoring_diff_nuclide = scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO
case (SCORE_KAPPA_FISSION)
scoring_diff_nuclide = scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % kappa_fission &
/= ZERO
case (SCORE_KEFF)
scoring_diff_nuclide = scoring_diff_nuclide .and. &
micro_xs(t % deriv % diff_nuclide) % nu_fission /= ZERO
@ -874,24 +861,6 @@ contains
end select
case (SCORE_KAPPA_FISSION)
select case (t % deriv % dep_var)
case (DIFF_NUCLIDE_DENSITY)
if (i_nuclide == -1 .and. &
materials(p % material)%id== t % deriv % diff_material) then
score = score * (t % deriv % accumulator &
+ micro_xs(t % deriv % diff_nuclide) % kappa_fission &
/ material_xs % kappa_fission)
else if (scoring_diff_nuclide) then
score = score * (t % deriv % accumulator + ONE &
/ atom_density)
else
score = score * t % deriv % accumulator
end if
end select
case (SCORE_KEFF)
select case (t % deriv % dep_var)
@ -901,7 +870,7 @@ contains
case (DIFF_NUCLIDE_DENSITY)
if (scoring_diff_nuclide) then
score = score * (t % deriv % accumulator + ONE &
/ mat % atom_density(l))
/ materials(p % material) % atom_density(l))
else
score = score * t % deriv % accumulator
end if
@ -1185,9 +1154,8 @@ contains
!===============================================================================
subroutine score_fission_eout(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), pointer :: t
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
integer :: i ! index of outgoing energy filter
@ -1518,9 +1486,9 @@ contains
logical :: end_in_mesh ! ending coordinates inside mesh?
real(8) :: theta
real(8) :: phi
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(Material), pointer :: mat
type(Material), pointer :: mat
t => tallies(i_tally)
matching_bins(1:t%n_filters) = 1
@ -1894,7 +1862,7 @@ contains
integer :: offset ! offset for distribcell
real(8) :: E ! particle energy
real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
found_bin = .true.
@ -2118,7 +2086,7 @@ contains
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
TALLY_LOOP: do i = 1, active_current_tallies % size()

82
src/tally_header.F90
View file

@ -58,10 +58,6 @@ module tally_header
integer :: offset = 0 ! Only used for distribcell filters
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
! Type-Bound procedures
contains
procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter
end type TallyFilter
@ -144,84 +140,6 @@ module tally_header
! Derivative for differentially tallies
type(TallyDerivative), allocatable :: deriv
! Type-Bound procedures
contains
procedure :: clear => tallyobject_clear ! Deallocates TallyObject
end type TallyObject
contains
!===============================================================================
! TALLYFILTER_CLEAR deallocates a TallyFilter element and sets it to its as
! initialized state.
!===============================================================================
subroutine tallyfilter_clear(this)
class(TallyFilter), intent(inout) :: this ! The TallyFilter to be cleared
this % type = NONE
this % n_bins = 0
if (allocated(this % int_bins)) &
deallocate(this % int_bins)
if (allocated(this % real_bins)) &
deallocate(this % real_bins)
end subroutine tallyfilter_clear
!===============================================================================
! TALLYOBJECT_CLEAR deallocates a TallyObject element and sets it to its as
! initialized state.
!===============================================================================
subroutine tallyobject_clear(this)
class(TallyObject), intent(inout) :: this ! The TallyObject to be cleared
integer :: i ! Loop Index
! This routine will go through each item in TallyObject and set the value
! to its default, as-initialized values, including deallocations.
this % name = ""
if (allocated(this % filters)) then
do i = 1, size(this % filters)
call this % filters(i) % clear()
end do
deallocate(this % filters)
end if
if (allocated(this % stride)) &
deallocate(this % stride)
this % find_filter = 0
this % n_nuclide_bins = 0
if (allocated(this % nuclide_bins)) &
deallocate(this % nuclide_bins)
this % all_nuclides = .false.
this % n_score_bins = 0
if (allocated(this % score_bins)) &
deallocate(this % score_bins)
if (allocated(this % moment_order)) &
deallocate(this % moment_order)
this % n_user_score_bins = 0
if (allocated(this % results)) &
deallocate(this % results)
this % reset = .false.
this % n_realizations = 0
if (allocated(this % triggers)) &
deallocate (this % triggers)
this % n_triggers = 0
if (allocated(this % deriv)) &
deallocate (this % deriv)
end subroutine tallyobject_clear
end module tally_header

4
src/tally_initialize.F90
View file

@ -38,7 +38,7 @@ contains
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: max_n_filters = 0 ! maximum number of filters
type(TallyObject), pointer :: t => null()
type(TallyObject), pointer :: t
TALLY_LOOP: do i = 1, n_tallies
! Get pointer to tally
@ -88,7 +88,7 @@ contains
integer :: k ! loop index for bins
integer :: bin ! filter bin entries
integer :: type ! type of tally filter
type(TallyObject), pointer :: t => null()
type(TallyObject), pointer :: t
! allocate tally map array -- note that we don't need a tally map for the
! energy_in and energy_out filters

12
src/timer_header.F90
View file

@ -29,13 +29,11 @@ contains
!===============================================================================
subroutine timer_start(self)
class(Timer), intent(inout) :: self
! Turn timer on and measure starting time
self % running = .true.
call system_clock(self % start_counts)
end subroutine timer_start
!===============================================================================
@ -43,7 +41,6 @@ contains
!===============================================================================
function timer_get_value(self) result(elapsed)
class(Timer), intent(in) :: self ! the timer
real(8) :: elapsed ! total elapsed time
@ -58,7 +55,6 @@ contains
else
elapsed = self % elapsed
end if
end function timer_get_value
!===============================================================================
@ -66,30 +62,26 @@ contains
!===============================================================================
subroutine timer_stop(self)
class(Timer), intent(inout) :: self
! Check to make sure timer was running
if (.not. self % running) return
! Stop timer and add time
self % elapsed = timer_get_value(self)
self % elapsed = self % get_value()
self % running = .false.
end subroutine timer_stop
!===============================================================================
! TIMER_RESET resets a timer to have a zero value
!===============================================================================
subroutine timer_reset(self)
pure subroutine timer_reset(self)
class(Timer), intent(inout) :: self
self % running = .false.
self % start_counts = 0
self % elapsed = ZERO
end subroutine timer_reset
end module timer_header

10
src/trigger_header.F90
View file

@ -1,6 +1,6 @@
module trigger_header
use constants, only: NONE, N_FILTER_TYPES
use constants, only: NONE, N_FILTER_TYPES, ZERO
implicit none
@ -13,9 +13,9 @@ module trigger_header
real(8) :: threshold ! a convergence threshold
character(len=52) :: score_name ! the name of the score
integer :: score_index ! the index of the score
real(8) :: variance=0.0 ! temp variance container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: rel_err =0.0 ! temp rel. err. container
real(8) :: variance = ZERO ! temp variance container
real(8) :: std_dev = ZERO ! temp std. dev. container
real(8) :: rel_err = ZERO ! temp rel. err. container
end type TriggerObject
!===============================================================================
@ -23,7 +23,7 @@ module trigger_header
!===============================================================================
type KTrigger
integer :: trigger_type = 0
real(8) :: threshold = 0
real(8) :: threshold = ZERO
end type KTrigger
end module trigger_header

2
src/xml_interface.F90
View file

@ -117,7 +117,7 @@ contains
type(Node), pointer, intent(out) :: out_ptr
logical :: found_
type(NodeList), pointer :: elem_list => null()
type(NodeList), pointer :: elem_list
! Set found to false
found_ = .false.

5
tests/test_complex_cell/geometry.xml
View file

@ -15,10 +15,11 @@
<surface id="14" type="y-plane" coeffs="4" />
<surface id="15" type="y-plane" coeffs="7" />
<surface id="16" type="y-plane" coeffs="10" boundary="vacuum" />
<surface id="17" type="y-plane" coeffs="0" />
<cell id="1" material="1" region="3 -4 13 -14" />
<cell id="2" material="2" region="2 -5 12 -15 ~(3 -4 13 -14)" />
<cell id="3" material="3" region="7 -6 11 -16 (-2 | 5 | -12 | 15)" />
<cell id="4" material="4" region="((1 -7) 11 -16) ~(2 -5 (12 -15))" />
<cell id="3" material="3" region="((1 -7 17 -16) | (7 -6 11 -17)) (-2 | 5 | -12 | 15)" />
<cell id="4" material="4" region="((1 -7 11 -17) | (7 -6 17 -16)) ~(2 -5 12 -15)" />
</geometry>

18
tests/test_complex_cell/results_true.dat
View file

@ -1,11 +1,11 @@
k-combined:
2.651570E-01 2.116381E-03
2.565769E-01 8.980879E-04
tally 1:
2.639097E+00
1.394398E+00
2.743740E+00
1.506124E+00
1.041248E+00
2.177204E-01
1.087210E-01
2.365126E-03
2.584080E+00
1.335682E+00
2.763580E+00
1.528633E+00
1.007148E+00
2.031543E-01
1.113696E-01
2.485351E-03

4
tests/test_many_scores/results_true.dat
View file

@ -33,8 +33,8 @@ tally 1:
7.620560E-01
1.816851E+00
1.102658E+00
1.338067E+02
5.986137E+03
1.337996E+02
5.985519E+03
2.247257E+01
1.683779E+02
1.512960E-01

20
tests/test_score_kappafission/results_true.dat
View file

@ -1,17 +1,17 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.266048E+02
1.049833E+04
2.266169E+02
1.049923E+04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.366139E+02
3.859561E+03
1.366590E+02
3.861651E+03
tally 2:
2.402814E+02
1.174003E+04
2.403775E+02
1.175130E+04
0.000000E+00
0.000000E+00
0.000000E+00
@ -19,11 +19,11 @@ tally 2:
1.270420E+02
3.297537E+03
tally 3:
2.217075E+02
1.003168E+04
2.217588E+02
1.003581E+04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.375693E+02
3.872389E+03
1.376303E+02
3.875598E+03