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

Merge pull request #4 from paulromano/multitemp-refactor

Browse source 
Multi-temperature data refactor
This commit is contained in:
Adam Nelson 2016-09-02 14:50:38 -04:00 • committed by GitHub
commit b56ceaa0f6
157 changed files with 4251 additions and 3885 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
docs/source/_images/openmc.png
View file

Binary file not shown.
Side by side

Before

Width:  |  Height:  |  Size: 8.7 KiB

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

Binary file not shown.
Side by side

Before

Width:  |  Height:  |  Size: 6.8 KiB

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

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 15 KiB

60
docs/source/_images/openmc_logo.svg Normal file
View file

@ -0,0 +1,60 @@
<?xml version="1.0" encoding="utf-8"?>
<!-- Generator: Adobe Illustrator 16.0.0, SVG Export Plug-In . SVG Version: 6.00 Build 0) -->
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg version="1.1" id="Layer_1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
width="257.157px" height="60px" viewBox="0 0 257.157 60" enable-background="new 0 0 257.157 60" xml:space="preserve">
<g>
<g>
<path fill="#656565" d="M81.676,48.809c-2.67,0-5.128-0.469-7.367-1.409c-2.243-0.937-4.179-2.206-5.81-3.805
c-1.632-1.601-2.901-3.477-3.807-5.638c-0.908-2.158-1.36-4.475-1.36-6.947v-0.099c0-2.472,0.46-4.788,1.385-6.945
c0.922-2.159,2.199-4.057,3.832-5.688c1.63-1.631,3.576-2.916,5.834-3.856s4.724-1.409,7.392-1.409
c2.67,0,5.125,0.469,7.367,1.409s4.179,2.21,5.812,3.808c1.63,1.6,2.899,3.479,3.805,5.637c0.907,2.161,1.361,4.476,1.361,6.947
v0.099c0,2.472-0.463,4.79-1.385,6.945c-0.922,2.161-2.2,4.059-3.832,5.688c-1.63,1.629-3.576,2.918-5.834,3.854
C86.81,48.34,84.344,48.809,81.676,48.809z M81.773,41.788c1.517,0,2.918-0.279,4.203-0.839c1.286-0.561,2.382-1.336,3.288-2.325
c0.907-0.986,1.614-2.131,2.126-3.437c0.512-1.302,0.768-2.694,0.768-4.178v-0.099c0-1.481-0.256-2.883-0.768-4.202
c-0.512-1.318-1.235-2.475-2.175-3.461c-0.94-0.991-2.053-1.771-3.338-2.35c-1.285-0.576-2.687-0.864-4.201-0.864
c-1.552,0-2.958,0.279-4.228,0.838c-1.27,0.562-2.358,1.336-3.264,2.327c-0.907,0.987-1.616,2.132-2.125,3.436
c-0.511,1.303-0.768,2.693-0.768,4.178v0.099c0,1.483,0.257,2.884,0.768,4.205c0.51,1.318,1.234,2.47,2.175,3.46
c0.939,0.988,2.042,1.772,3.313,2.347C78.815,41.499,80.224,41.788,81.773,41.788z"/>
<path fill="#656565" d="M102.145,21.715h7.516v3.808c0.924-1.253,2.035-2.283,3.338-3.09c1.301-0.809,2.942-1.211,4.92-1.211
c1.549,0,3.048,0.297,4.498,0.889c1.453,0.594,2.737,1.475,3.859,2.645c1.119,1.17,2.019,2.604,2.694,4.301
c0.675,1.7,1.013,3.652,1.013,5.859v0.1c0,2.209-0.338,4.162-1.013,5.858c-0.675,1.697-1.565,3.133-2.67,4.301
c-1.105,1.173-2.381,2.054-3.832,2.648c-1.452,0.594-2.966,0.889-4.548,0.889c-2.012,0-3.667-0.397-4.971-1.187
c-1.302-0.79-2.396-1.712-3.287-2.768v11.371h-7.516V21.715z M115.989,42.332c0.889,0,1.722-0.172,2.498-0.518
c0.773-0.348,1.458-0.843,2.051-1.484c0.592-0.641,1.064-1.408,1.409-2.298c0.347-0.89,0.52-1.896,0.52-3.018v-0.1
c0-1.087-0.173-2.082-0.52-2.99c-0.345-0.907-0.816-1.682-1.409-2.325c-0.593-0.642-1.278-1.137-2.051-1.482
c-0.776-0.345-1.609-0.52-2.498-0.52s-1.722,0.175-2.496,0.52c-0.775,0.346-1.451,0.841-2.028,1.482
c-0.577,0.644-1.037,1.418-1.385,2.325c-0.345,0.908-0.518,1.903-0.518,2.99v0.1c0,1.087,0.173,2.086,0.518,2.991
c0.348,0.906,0.808,1.684,1.385,2.324c0.577,0.642,1.253,1.137,2.028,1.484C114.267,42.16,115.101,42.332,115.989,42.332z"/>
<path fill="#656565" d="M145.62,48.809c-1.979,0-3.814-0.329-5.514-0.986c-1.696-0.661-3.162-1.6-4.399-2.816
c-1.237-1.223-2.199-2.666-2.891-4.331c-0.693-1.661-1.041-3.517-1.041-5.559v-0.102c0-1.877,0.321-3.658,0.964-5.34
c0.645-1.682,1.543-3.147,2.695-4.4c1.154-1.251,2.53-2.241,4.129-2.967c1.6-0.725,3.37-1.086,5.316-1.086
c2.206,0,4.118,0.394,5.733,1.187c1.617,0.791,2.958,1.854,4.03,3.188c1.072,1.335,1.864,2.867,2.374,4.598
c0.511,1.731,0.766,3.537,0.766,5.417c0,0.295-0.008,0.61-0.024,0.938c-0.016,0.331-0.04,0.676-0.074,1.038h-18.442
c0.364,1.716,1.112,3.009,2.25,3.881c1.138,0.873,2.547,1.313,4.228,1.313c1.251,0,2.373-0.216,3.363-0.646
c0.988-0.427,2.01-1.118,3.064-2.076l4.304,3.81c-1.254,1.547-2.771,2.76-4.552,3.633C150.12,48.372,148.026,48.809,145.62,48.809
z M150.464,32.89c-0.227-1.681-0.821-3.042-1.777-4.08c-0.956-1.037-2.226-1.557-3.807-1.557c-1.582,0-2.861,0.512-3.832,1.53
c-0.973,1.023-1.607,2.393-1.904,4.106H150.464z"/>
<path fill="#656565" d="M159.415,21.715h7.518v3.787c0.428-0.562,0.896-1.103,1.408-1.616c0.512-0.517,1.08-0.971,1.706-1.371
c0.625-0.397,1.318-0.712,2.078-0.946c0.757-0.233,1.613-0.347,2.569-0.347c2.867,0,5.084,0.873,6.65,2.619
c1.565,1.746,2.35,4.152,2.35,7.219v17.156h-7.518V33.471c0-1.776-0.395-3.118-1.186-4.021c-0.792-0.903-1.913-1.355-3.361-1.355
c-1.451,0-2.598,0.452-3.437,1.355c-0.841,0.903-1.261,2.245-1.261,4.021v14.745h-7.518V21.715z"/>
<path fill="#A31F34" d="M187.056,13.607h8.209l9.095,14.634l9.099-14.634h8.207v34.608h-7.513V25.619l-9.742,14.784h-0.196
l-9.645-14.634v22.446h-7.514V13.607z"/>
<path fill="#A31F34" d="M242.235,48.809c-2.537,0-4.893-0.459-7.071-1.385c-2.175-0.922-4.052-2.18-5.636-3.78
c-1.583-1.601-2.819-3.485-3.709-5.66c-0.888-2.179-1.333-4.501-1.333-6.974v-0.099c0-2.472,0.445-4.788,1.333-6.945
c0.89-2.159,2.126-4.057,3.709-5.688c1.584-1.631,3.477-2.916,5.687-3.856c2.207-0.94,4.647-1.409,7.317-1.409
c1.613,0,3.091,0.132,4.425,0.396c1.336,0.265,2.547,0.625,3.635,1.089c1.088,0.461,2.093,1.021,3.016,1.682
c0.923,0.657,1.781,1.385,2.57,2.175l-4.846,5.586c-1.349-1.218-2.728-2.175-4.128-2.866c-1.401-0.692-2.974-1.038-4.721-1.038
c-1.453,0-2.793,0.279-4.03,0.838c-1.235,0.562-2.299,1.336-3.188,2.327c-0.89,0.987-1.582,2.132-2.077,3.436
c-0.494,1.303-0.742,2.693-0.742,4.178v0.099c0,1.483,0.248,2.884,0.742,4.205c0.495,1.318,1.178,2.47,2.053,3.46
c0.872,0.988,1.929,1.772,3.164,2.347c1.235,0.576,2.594,0.865,4.079,0.865c1.979,0,3.65-0.362,5.02-1.086
c1.366-0.726,2.726-1.714,4.079-2.968l4.844,4.895c-0.89,0.958-1.814,1.814-2.768,2.571c-0.958,0.758-2.003,1.411-3.142,1.953
c-1.136,0.544-2.379,0.959-3.733,1.236C245.433,48.669,243.915,48.809,242.235,48.809z"/>
</g>
<path fill="#A31F34" d="M30.614,1.024c-10.866,0-20.325,5.991-25.284,14.845h20.495l13.167,22.804l-11.609,20.11
c1.062,0.119,2.139,0.192,3.231,0.192c16.003,0,28.975-12.972,28.975-28.976C59.589,13.997,46.617,1.024,30.614,1.024z"/>
<path fill="#656565" d="M16.416,55.243h5.807l9.567-16.57l-9.567-16.566H3.093l-0.607,1.052C1.954,25.355,1.639,27.638,1.639,30
C1.639,40.84,7.603,50.274,16.416,55.243z"/>
</g>
</svg>
Side by side

After

Width:  |  Height:  |  Size: 5.9 KiB

4
docs/source/_static/theme_overrides.css
View file

@ -16,3 +16,7 @@
.wy-table, .rst-content table.docutils, .rst-content table.field-list {
margin-bottom: 0px;
}
.wy-side-nav-search {
background-color: #343131;
}

2
docs/source/conf.py
View file

@ -129,7 +129,7 @@ if not on_rtd:
html_theme = 'sphinx_rtd_theme'
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
html_logo = '_images/openmc200px.png'
html_logo = '_images/openmc_logo.png'
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".

23
docs/source/io_formats/nuclear_data.rst
View file

@ -26,8 +26,9 @@ Incident Neutron Data
temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Datasets: - **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
**/<nuclide name>/reactions/reaction_<mt>/**
@ -44,11 +45,14 @@ temperature-dependent data set. For example, the data set corresponding to
temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Attributes: - **threshold_idx** (*int*) -- Index on the energy grid that the
reaction threshold corresponds to for temperature TTT (in Kelvin)
:Datasets:
- **xs** (*double[]*) -- Cross section values tabulated against the
nuclide energy grid for temperature TTT (in Kelvin)
:Datasets: - **xs** (*double[]*) -- Cross section values tabulated against the
nuclide energy grid for temperature TTT (in Kelvin)
:Attributes:
- **threshold_idx** (*int*) -- Index on the energy
grid that the reaction threshold corresponds to for
temperature TTT (in Kelvin)
**/<nuclide name>/reactions/reaction_<mt>/product_<j>/**
@ -112,7 +116,7 @@ Thermal Neutron Scattering Data
**/<thermal name>/**
:Attributes: - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses
- **zaids** (*int[]*) -- ZAID identifiers for which the thermal
- **nuclides** (*char[][]*) -- Names of nuclides for which the thermal
scattering data applies to
- **secondary_mode** (*char[]*) -- Indicates how the inelastic
outgoing angle-energy distributions are represented ('equal',
@ -124,8 +128,9 @@ Thermal Neutron Scattering Data
temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Datasets: - **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin)
**/<thermal name>/elastic/<TTT>K/**

37
docs/source/methods/cross_sections.rst
View file

@ -53,12 +53,12 @@ speed up the calculation.
Logarithmic Mapping
+++++++++++++++++++
To speed up energy grid searches, OpenMC uses logarithmic mapping technique
[Brown]_ to limit the range of energies that must be searched for each
nuclide. The entire energy range is divided up into equal-lethargy segments, and
the bounding energies of each segment are mapped to bounding indices on each of
the nuclide energy grids. By default, OpenMC uses 8000 equal-lethargy segments
as recommended by Brown.
To speed up energy grid searches, OpenMC uses a `logarithmic mapping technique`_
to limit the range of energies that must be searched for each nuclide. The
entire energy range is divided up into equal-lethargy segments, and the bounding
energies of each segment are mapped to bounding indices on each of the nuclide
energy grids. By default, OpenMC uses 8000 equal-lethargy segments as
recommended by Brown.
Other Methods
+++++++++++++
@ -74,9 +74,9 @@ offers support for an experimental data format called windowed multipole (WMP).
This data format requires less memory than pointwise cross sections, and it
allows on-the-fly Doppler broadening to arbitrary temperature.
The multipole method was introduced by [Hwang]_ and the faster windowed
multipole method by [Josey]_. In the multipole format, cross section resonances
are represented by poles, :math:`p_j`, and residues, :math:`r_j`, in the complex
The multipole method was introduced by Hwang_ and the faster windowed multipole
method by Josey_. In the multipole format, cross section resonances are
represented by poles, :math:`p_j`, and residues, :math:`r_j`, in the complex
plane. The 0K cross sections in the resolved resonance region can be computed
by summing up a contribution from each pole:
@ -232,21 +232,10 @@ sections. This allows flexibility for the model to use highly anisotropic
scattering information in the water while the fuel can be simulated with linear
or even isotropic scattering.
.. only:: html
.. rubric:: References
.. [Brown] Forrest B. Brown, "New Hash-based Energy Lookup Algorithm for Monte
Carlo codes," LA-UR-14-24530, Los Alamos National Laboratory (2014).
.. [Hwang] R. N. Hwang, "A Rigorous Pole Representation of Multilevel Cross
Sections and Its Practical Application," *Nucl. Sci. Eng.*, **96**,
192-209 (1987).
.. [Josey] Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed
Multipole for Cross Section Doppler Broadening," *J. Comp. Phys*,
**307**, 715-727 (2016). http://dx.doi.org/10.1016/j.jcp.2015.08.013
.. _logarithmic mapping technique:
https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
.. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381
.. _Josey: http://dx.doi.org/10.1016/j.jcp.2015.08.013
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi
.. _NJOY: http://t2.lanl.gov/codes.shtml

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

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

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

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

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

@ -214,7 +214,6 @@
"source": [
"# Instantiate a Materials collection and export to XML\n",
"materials_file = openmc.Materials([inf_medium])\n",
"materials_file.default_xs = '71c'\n",
"materials_file.export_to_xml()"
]
},
@ -499,23 +498,37 @@
"output_type": "stream",
"text": [
"\n",
" .d88888b. 888b d888 .d8888b.\n",
" d88P\" \"Y88b 8888b d8888 d88P Y88b\n",
" 888 888 88888b.d88888 888 888\n",
" 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 \n",
" 888 888 888 \"88b d8P Y8b 888 \"88b 888 Y888P 888 888 \n",
" 888 888 888 888 88888888 888 888 888 Y8P 888 888 888\n",
" Y88b. .d88P 888 d88P Y8b. 888 888 888 \" 888 Y88b d88P\n",
" \"Y88888P\" 88888P\" \"Y8888 888 888 888 888 \"Y8888P\"\n",
"__________________888______________________________________________________\n",
" 888\n",
" 888\n",
" %%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ################## %%%%%%%%%%%%%%%%%%%%%%%\n",
" ################### %%%%%%%%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%%%%%%\n",
" ##################### %%%%%%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%\n",
" ################# %%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%\n",
" ############ %%%%%%%%%%%%%%%\n",
" ######## %%%%%%%%%%%%%%\n",
" %%%%%%%%%%%\n",
"\n",
" Copyright: 2011-2016 Massachusetts Institute of Technology\n",
" License: http://openmc.readthedocs.io/en/latest/license.html\n",
" Version: 0.7.1\n",
" Git SHA1: 3d68c07625e33cd64188df03ee03e9c31b3d4b74\n",
" Date/Time: 2016-07-22 21:03:18\n",
" | The OpenMC Monte Carlo Code\n",
" Copyright | 2011-2016 Massachusetts Institute of Technology\n",
" License | http://openmc.readthedocs.io/en/latest/license.html\n",
" Version | 0.8.0\n",
" Git SHA1 | fbebf7bf709fe2fe1813af95bff9b29c0d59312c\n",
" Date/Time | 2016-08-31 10:40:13\n",
" OpenMP Threads | 4\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
@ -525,12 +538,12 @@
" Reading geometry XML file...\n",
" Reading cross sections XML file...\n",
" Reading materials XML file...\n",
" Reading H1.71c from /home/romano/openmc/data/nndc_hdf5/H1_71c.h5\n",
" Reading O16.71c from /home/romano/openmc/data/nndc_hdf5/O16_71c.h5\n",
" Reading U235.71c from /home/romano/openmc/data/nndc_hdf5/U235_71c.h5\n",
" Reading U238.71c from /home/romano/openmc/data/nndc_hdf5/U238_71c.h5\n",
" Reading Zr90.71c from /home/romano/openmc/data/nndc_hdf5/Zr90_71c.h5\n",
" Maximum neutron transport energy: 20.0000 MeV for H1.71c\n",
" Reading H1 from /home/romano/openmc/data/nndc_hdf5/H1.h5\n",
" Reading O16 from /home/romano/openmc/data/nndc_hdf5/O16.h5\n",
" Reading U235 from /home/romano/openmc/data/nndc_hdf5/U235.h5\n",
" Reading U238 from /home/romano/openmc/data/nndc_hdf5/U238.h5\n",
" Reading Zr90 from /home/romano/openmc/data/nndc_hdf5/Zr90.h5\n",
" Maximum neutron transport energy: 20.0000 MeV for H1\n",
" Reading tallies XML file...\n",
" Building neighboring cells lists for each surface...\n",
" Initializing source particles...\n",
@ -600,20 +613,20 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 3.2300E-01 seconds\n",
" Reading cross sections = 1.6900E-01 seconds\n",
" Total time in simulation = 1.9882E+01 seconds\n",
" Time in transport only = 1.9869E+01 seconds\n",
" Time in inactive batches = 2.6590E+00 seconds\n",
" Time in active batches = 1.7223E+01 seconds\n",
" Total time for initialization = 3.9900E-01 seconds\n",
" Reading cross sections = 2.6500E-01 seconds\n",
" Total time in simulation = 1.1488E+01 seconds\n",
" Time in transport only = 1.1152E+01 seconds\n",
" Time in inactive batches = 1.2180E+00 seconds\n",
" Time in active batches = 1.0270E+01 seconds\n",
" Time synchronizing fission bank = 4.0000E-03 seconds\n",
" Sampling source sites = 4.0000E-03 seconds\n",
" SEND/RECV source sites = 0.0000E+00 seconds\n",
" Sampling source sites = 3.0000E-03 seconds\n",
" SEND/RECV source sites = 1.0000E-03 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 0.0000E+00 seconds\n",
" Total time elapsed = 2.0217E+01 seconds\n",
" Calculation Rate (inactive) = 9402.03 neutrons/second\n",
" Calculation Rate (active) = 5806.19 neutrons/second\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 1.1901E+01 seconds\n",
" Calculation Rate (inactive) = 20525.5 neutrons/second\n",
" Calculation Rate (active) = 9737.10 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
@ -894,7 +907,7 @@
" <td>6.250000e-07</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>-3.774758e-15</td>\n",
" <td>-2.886580e-15</td>\n",
" <td>0.011292</td>\n",
" </tr>\n",
" <tr>\n",
@ -904,7 +917,7 @@
" <td>2.000000e+01</td>\n",
" <td>total</td>\n",
" <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
" <td>1.443290e-15</td>\n",
" <td>-5.551115e-16</td>\n",
" <td>0.002570</td>\n",
" </tr>\n",
" </tbody>\n",
@ -917,8 +930,8 @@
"1 1 6.25e-07 2.00e+01 total \n",
"\n",
" score mean std. dev. \n",
"0 (((total / flux) - (absorption / flux)) - (sca... -3.77e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... 1.44e-15 2.57e-03 "
"0 (((total / flux) - (absorption / flux)) - (sca... -2.89e-15 1.13e-02 \n",
"1 (((total / flux) - (absorption / flux)) - (sca... -5.55e-16 2.57e-03 "
]
},
"execution_count": 22,
@ -1167,21 +1180,21 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"display_name": "Python 3",
"language": "python",
"name": "python2"
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

316
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

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

File diff suppressed because one or more lines are too long

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

@ -105,7 +105,6 @@
"source": [
"# Instantiate a Materials collection\n",
"materials_file = openmc.Materials((fuel, water, zircaloy))\n",
"materials_file.default_xs = '71c'\n",
"\n",
"# Export to \"materials.xml\"\n",
"materials_file.export_to_xml()"
@ -339,7 +338,7 @@
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AHFwInLqDpadAAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDctMjJUMjE6Mzk6\nNDYtMDU6MDBOOEOsAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA3LTIyVDIxOjM5OjQ2LTA1OjAw\nP2X7EAAAAABJRU5ErkJggg==\n",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEX///9yEhLpgJFNv8Tq\nQYT7AAAAAWJLR0QAiAUdSAAAAAd0SU1FB+AIHw8dMt59x4sAAALKSURBVGje7dpLcqQwDAbgHHE2\nYeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmNP+HDhw8fPnz48Kf6VH9G\n+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4zPji99z0/AJ4n1lfvJ6f\nnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6pA0wfln+ho/fwgYYn19C\n/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tNDbSGz7T0SBEWw4vLXzbQ\n6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X58wZaxWd1+fMGiuFvir8b\nvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV873hB8UnM3xzANtf8nb4\ndwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7T/ppARBvp48UwJnelT5S\nACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4//Jve+fhsH6Ctv7n8PTzj\nvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V32/o9+fl389Xnx+g5x/o\n+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6/4Le/6D3T/D9V67Y/ZsV\nQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/gPs/0P4TtP8F7r9J3AIO\n9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTuf4X7b+H+X7T/+BPuf3aM\n8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIwMTYtMDgtMzFUMTA6Mjk6\nNTAtMDU6MDBsyrzpAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE2LTA4LTMxVDEwOjI5OjUwLTA1OjAw\nHZcEVQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
@ -527,23 +526,37 @@
"output_type": "stream",
"text": [
"\n",
" .d88888b. 888b d888 .d8888b.\n",
" d88P\" \"Y88b 8888b d8888 d88P Y88b\n",
" 888 888 88888b.d88888 888 888\n",
" 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 \n",
" 888 888 888 \"88b d8P Y8b 888 \"88b 888 Y888P 888 888 \n",
" 888 888 888 888 88888888 888 888 888 Y8P 888 888 888\n",
" Y88b. .d88P 888 d88P Y8b. 888 888 888 \" 888 Y88b d88P\n",
" \"Y88888P\" 88888P\" \"Y8888 888 888 888 888 \"Y8888P\"\n",
"__________________888______________________________________________________\n",
" 888\n",
" 888\n",
" %%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%%%%%%%%%\n",
" ################## %%%%%%%%%%%%%%%%%%%%%%%\n",
" ################### %%%%%%%%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%%%%%%\n",
" ##################### %%%%%%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%%\n",
" ####################### %%%%%%%%%%%%%%%%%\n",
" ###################### %%%%%%%%%%%%%%%%%\n",
" #################### %%%%%%%%%%%%%%%%%\n",
" ################# %%%%%%%%%%%%%%%%%\n",
" ############### %%%%%%%%%%%%%%%%\n",
" ############ %%%%%%%%%%%%%%%\n",
" ######## %%%%%%%%%%%%%%\n",
" %%%%%%%%%%%\n",
"\n",
" Copyright: 2011-2016 Massachusetts Institute of Technology\n",
" License: http://openmc.readthedocs.io/en/latest/license.html\n",
" Version: 0.7.1\n",
" Git SHA1: 3d68c07625e33cd64188df03ee03e9c31b3d4b74\n",
" Date/Time: 2016-07-22 21:39:46\n",
" | The OpenMC Monte Carlo Code\n",
" Copyright | 2011-2016 Massachusetts Institute of Technology\n",
" License | http://openmc.readthedocs.io/en/latest/license.html\n",
" Version | 0.8.0\n",
" Git SHA1 | fbebf7bf709fe2fe1813af95bff9b29c0d59312c\n",
" Date/Time | 2016-08-31 10:29:51\n",
" OpenMP Threads | 4\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
@ -553,13 +566,13 @@
" Reading geometry XML file...\n",
" Reading cross sections XML file...\n",
" Reading materials XML file...\n",
" Reading U235.71c from /home/romano/openmc/data/nndc_hdf5/U235_71c.h5\n",
" Reading U238.71c from /home/romano/openmc/data/nndc_hdf5/U238_71c.h5\n",
" Reading O16.71c from /home/romano/openmc/data/nndc_hdf5/O16_71c.h5\n",
" Reading H1.71c from /home/romano/openmc/data/nndc_hdf5/H1_71c.h5\n",
" Reading B10.71c from /home/romano/openmc/data/nndc_hdf5/B10_71c.h5\n",
" Reading Zr90.71c from /home/romano/openmc/data/nndc_hdf5/Zr90_71c.h5\n",
" Maximum neutron transport energy: 20.0000 MeV for U235.71c\n",
" Reading U235 from /home/romano/openmc/data/nndc_hdf5/U235.h5\n",
" Reading U238 from /home/romano/openmc/data/nndc_hdf5/U238.h5\n",
" Reading O16 from /home/romano/openmc/data/nndc_hdf5/O16.h5\n",
" Reading H1 from /home/romano/openmc/data/nndc_hdf5/H1.h5\n",
" Reading B10 from /home/romano/openmc/data/nndc_hdf5/B10.h5\n",
" Reading Zr90 from /home/romano/openmc/data/nndc_hdf5/Zr90.h5\n",
" Maximum neutron transport energy: 20.0000 MeV for U235\n",
" Reading tallies XML file...\n",
" Building neighboring cells lists for each surface...\n",
" Initializing source particles...\n",
@ -599,20 +612,20 @@
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 3.5600E-01 seconds\n",
" Reading cross sections = 2.3400E-01 seconds\n",
" Total time in simulation = 1.8333E+01 seconds\n",
" Time in transport only = 1.8325E+01 seconds\n",
" Time in inactive batches = 2.6950E+00 seconds\n",
" Time in active batches = 1.5638E+01 seconds\n",
" Time synchronizing fission bank = 1.0000E-03 seconds\n",
" Sampling source sites = 0.0000E+00 seconds\n",
" SEND/RECV source sites = 1.0000E-03 seconds\n",
" Total time for initialization = 4.3100E-01 seconds\n",
" Reading cross sections = 3.0500E-01 seconds\n",
" Total time in simulation = 8.9870E+00 seconds\n",
" Time in transport only = 8.9500E+00 seconds\n",
" Time in inactive batches = 1.1950E+00 seconds\n",
" Time in active batches = 7.7920E+00 seconds\n",
" Time synchronizing fission bank = 5.0000E-03 seconds\n",
" Sampling source sites = 5.0000E-03 seconds\n",
" SEND/RECV source sites = 0.0000E+00 seconds\n",
" Time accumulating tallies = 0.0000E+00 seconds\n",
" Total time for finalization = 1.0000E-03 seconds\n",
" Total time elapsed = 1.8711E+01 seconds\n",
" Calculation Rate (inactive) = 4638.22 neutrons/second\n",
" Calculation Rate (active) = 2398.00 neutrons/second\n",
" Total time elapsed = 9.4370E+00 seconds\n",
" Calculation Rate (inactive) = 10460.3 neutrons/second\n",
" Calculation Rate (active) = 4812.63 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",

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

@ -334,6 +334,7 @@ Functions
:nosignatures:
openmc.model.create_triso_lattice
openmc.model.pack_trisos
--------------------------------------------
:mod:`openmc.data` -- Nuclear Data Interface

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

@ -281,6 +281,8 @@ based on the recommended value in LA-UR-14-24530_.
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
.. _multipole_library:
``<multipole_library>`` Element
-------------------------------
@ -290,8 +292,8 @@ OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range
cross sections. If this element is absent from the settings.xml file, the
:envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used.
.. note:: The <use_windowed_multipole> element must also be set to "true"
for windowed multipole functionality.
.. note:: The :ref:`temperature_method` must also be set to "multipole" for
windowed multipole functionality.
``<max_order>`` Element
---------------------------
@ -707,6 +709,47 @@ survival biasing, otherwise known as implicit capture or absorption.
*Default*: false
.. _temperature_default:
``<temperature_default>`` Element
---------------------------------
The ``<temperature_default>`` element specifies a default temperature in Kelvin
that is to be applied to cells in the absence of an explicit cell temperature or
a material default temperature.
*Default*: 293.6 K
.. _temperature_method:
``<temperature_method>`` Element
--------------------------------
The ``<temperature_method>`` element has an accepted value of "nearest" or
"interpolation". A value of "nearest" indicates that for each cell, the nearest
temperature at which cross sections are given is to be applied, within a given
tolerance (see :ref:`temperature_tolerance`). A value of "interpolation"
indicates that cross sections are to be interpolated between temperatures at
which nuclear data are present. A value of "multipole" indicates that the
windowed multipole method should be used to evaluate temperature-dependent cross
sections in the resolved resonance range (a :ref:`windowed multipole library
<multipole_library>` must also be available).
*Default*: "nearest"
.. _temperature_tolerance:
``<temperature_tolerance>`` Element
-----------------------------------
The ``<temperature_tolerance>`` element specifies a tolerance in Kelvin that is
to be applied when the "nearest" temperature method is used. For example, if a
cell temperature is 340 K and the tolerance is 15 K, then the closest
temperature in the range of 325 K to 355 K will be used to evaluate cross
sections.
*Default*: 10 K
``<threads>`` Element
---------------------
@ -1083,7 +1126,9 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
specified for the "distributed temperature" feature. This will give each
unique instance of the cell its own temperature.
*Default*: The temperature of the coldest nuclide in the cell's material(s)
*Default*: If a material default temperature is supplied, it is used. In the
absence of a material default temperature, the :ref:`global default
temperature <temperature_default>` is used.
:rotation:
If the cell is filled with a universe, this element specifies the angles in
@ -1289,11 +1334,12 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: ""
:temperature:
An element with no attributes which is used to set the default temperature
of the material in Kelvin.
An element with no attributes which is used to set the temperature of the
material. This element accepts a maximum 6-character string that indicates
the default temperature rounded to the nearest integer in units of Kelvin,
e.g. "294K".
*Default*: If a material default temperature is not given and a cell
temperature is not specified, the :ref:`global default temperature
<temperature_default>` is used.
:density:
An element with attributes/sub-elements called ``value`` and ``units``. The
@ -1409,19 +1455,6 @@ Each ``material`` element can have the following attributes or sub-elements:
.. _IUPAC Isotopic Compositions of the Elements 2009:
http://pac.iupac.org/publications/pac/pdf/2011/pdf/8302x0397.pdf
``<default_temperature>`` Element
------------------------
In some circumstances, the temperature may be the same for many or
all materials in a given problem. In this case, rather than specifying the
``<temperature>`` element on every material, a ``<default_temperature>``
element can be used to set the default material temperature for any material
without an explicitly provided temperature. This element has no attributes and
accepts a maximum 6-character string that indicates the default temperature
rounded to the nearest integer in units of Kelvin, e.g. "294K".
*Default*: None
------------------------------------
Tallies Specification -- tallies.xml
------------------------------------
@ -2102,8 +2135,8 @@ attributes or sub-elements. These are not used in "voxel" plots:
*Default*: None
:meshlines:
The ``meshlines`` sub-element allows for plotting the boundaries of
a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per
The ``meshlines`` sub-element allows for plotting the boundaries of a
regular mesh on top of a plot. Only one ``meshlines`` element is allowed per
``plot`` element, and it must contain as attributes or sub-elements a mesh
type and a linewidth. Optionally, a color may be specified for the overlay:

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

@ -33,7 +33,6 @@ fuel.add_nuclide(u235, 1.)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([moderator, fuel])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -38,7 +38,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([fuel1, fuel2, moderator])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -37,7 +37,6 @@ iron.add_nuclide(fe56, 1.)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([moderator, fuel, iron])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -32,7 +32,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials((moderator, fuel))
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -32,7 +32,6 @@ moderator.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([moderator, fuel])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -102,7 +102,6 @@ borated_water.add_s_alpha_beta('c_H_in_H2O')
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([uo2, helium, zircaloy, borated_water])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

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

@ -19,7 +19,7 @@ groups = openmc.mgxs.EnergyGroups(group_edges=[1E-11, 0.0635E-6, 10.0E-6,
1.0E-4, 1.0E-3, 0.5, 1.0, 20.0])
# Instantiate the 7-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2.300K', groups)
uo2_xsdata = openmc.XSdata('UO2', groups)
uo2_xsdata.order = 0
uo2_xsdata.total = [0.1779492, 0.3298048, 0.4803882, 0.5543674,
0.3118013, 0.3951678, 0.5644058]
@ -41,7 +41,7 @@ uo2_xsdata.nu_fission = [2.005998E-02, 2.027303E-03, 1.570599E-02,
uo2_xsdata.chi = [5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07,
0.0000E+00, 0.0000E+00, 0.0000E+00]
h2o_xsdata = openmc.XSdata('LWTR.300K', groups)
h2o_xsdata = openmc.XSdata('LWTR', groups)
h2o_xsdata.order = 0
h2o_xsdata.total = [0.15920605, 0.412969593, 0.59030986, 0.58435,
0.718, 1.2544497, 2.650379]
@ -66,8 +66,8 @@ mg_cross_sections_file.export_to_xml()
###############################################################################
# Instantiate some Macroscopic Data
uo2_data = openmc.Macroscopic('UO2', '300K')
h2o_data = openmc.Macroscopic('LWTR', '300K')
uo2_data = openmc.Macroscopic('UO2')
h2o_data = openmc.Macroscopic('LWTR')
# Instantiate some Materials and register the appropriate Macroscopic objects
uo2 = openmc.Material(material_id=1, name='UO2 fuel')
@ -80,7 +80,6 @@ water.add_macroscopic(h2o_data)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([uo2, water])
materials_file.default_xs = '300K'
materials_file.export_to_xml()

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

@ -25,7 +25,6 @@ fuel.add_nuclide(u235, 1.)
# Instantiate a Materials collection and export to XML
materials_file = openmc.Materials([fuel])
materials_file.default_temperature = '294K'
materials_file.export_to_xml()

2
examples/xml/basic/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

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

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
examples/xml/lattice/nested/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

2
examples/xml/lattice/simple/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

3
examples/xml/pincell/materials.xml
View file

@ -1,9 +1,6 @@
<?xml version="1.0"?>
<materials>
<!-- By default, use 294K cross sections -->
<default_temperature>294K</default_temperature>
<!--
Since O-18 is not present in ENDF/B-VII, it was necessary to combine the
atom densities for O-17 and O-18 in any materials containing Oxygen.

3
examples/xml/pincell_multigroup/materials.xml
View file

@ -1,8 +1,5 @@
<?xml version="1.0"?>
<materials>
<!-- Set default xs set to use 300K data -->
<default_xs>300K</default_xs>
<!-- UO2 -->
<material id="1">
<density units="macro" value="1.0" />

32
examples/xml/pincell_multigroup/mg_cross_sections.xml
View file

@ -11,8 +11,8 @@
-->
<xsdata>
<!-- Meta data for this data -->
<name>UO2.300K</name>
<alias>UO2.300K</alias>
<name>UO2</name>
<alias>UO2</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -67,8 +67,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX1.300K</name>
<alias>MOX1.300K</alias>
<name>MOX1</name>
<alias>MOX1</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -124,8 +124,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX2.300K</name>
<alias>MOX2.300K</alias>
<name>MOX2</name>
<alias>MOX2</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -180,8 +180,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>MOX3.300K</name>
<alias>MOX3.300K</alias>
<name>MOX3</name>
<alias>MOX3</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -236,8 +236,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>FC.300K</name>
<alias>FC.300K</alias>
<name>FC</name>
<alias>FC</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
@ -286,8 +286,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>GT.300K</name>
<alias>GT.300K</alias>
<name>GT</name>
<alias>GT</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
@ -318,8 +318,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>LWTR.300K</name>
<alias>LWTR.300K</alias>
<name>LWTR</name>
<alias>LWTR</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
@ -351,8 +351,8 @@
<xsdata>
<!-- Meta data for this data -->
<name>CR.300K</name>
<alias>CR.300K</alias>
<name>CR</name>
<alias>CR</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>

2
examples/xml/reflective/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

9
openmc/data/neutron.py
View file

@ -100,7 +100,7 @@ class IncidentNeutron(EqualityMixin):
Parameters
----------
name : str
Name of the table
Name of the nuclide using the GND naming convention
atomic_number : int
Number of protons in the nucleus
mass_number : int
@ -133,7 +133,7 @@ class IncidentNeutron(EqualityMixin):
metastable : int
Metastable state of the nucleus. A value of zero indicates ground state.
name : str
ZAID identifier of the table, e.g. 92235.
Name of the nuclide using the GND naming convention
reactions : collections.OrderedDict
Contains the cross sections, secondary angle and energy distributions,
and other associated data for each reaction. The keys are the MT values
@ -534,11 +534,6 @@ class IncidentNeutron(EqualityMixin):
# Assign temperature to the running list
kTs = [ace.temperature]
# If mass number hasn't been specified, make an educated guess
zaid, xs = ace.name.split('.')
name, element, Z, mass_number, metastable = \
_get_metadata(int(zaid), metastable_scheme)
data = cls(name, Z, mass_number, metastable,
ace.atomic_weight_ratio, kTs)

24
openmc/data/thermal.py
View file

@ -8,7 +8,7 @@ import h5py
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import K_BOLTZMANN
from .data import K_BOLTZMANN, ATOMIC_SYMBOL
from .ace import Table, get_table
from .angle_energy import AngleEnergy
from .function import Tabulated1D
@ -156,7 +156,7 @@ class ThermalScattering(EqualityMixin):
Parameters
----------
name : str
ZAID identifier of the table, e.g. lwtr.10t.
Name of the material using GND convention, e.g. c_H_in_H2O
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
kTs : Iterable of float
@ -174,7 +174,7 @@ class ThermalScattering(EqualityMixin):
Inelastic scattering cross section derived in the incoherent
approximation
name : str
Name of the table, e.g. lwtr.20t.
Name of the material using GND convention, e.g. c_H_in_H2O
temperatures : Iterable of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings of the temperature,
@ -182,8 +182,8 @@ class ThermalScattering(EqualityMixin):
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of MeV.
zaids : Iterable of int
ZAID identifiers that the thermal scattering data applies to
nuclides : Iterable of str
Nuclide names that the thermal scattering data applies to
"""
@ -200,7 +200,7 @@ class ThermalScattering(EqualityMixin):
self.inelastic_mu_out = {}
self.inelastic_dist = {}
self.secondary_mode = None
self.zaids = []
self.nuclides = []
def __repr__(self):
if hasattr(self, 'name'):
@ -226,7 +226,7 @@ class ThermalScattering(EqualityMixin):
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['zaids'] = self.zaids
g.attrs['nuclides'] = np.string_(self.nuclides)
g.attrs['secondary_mode'] = np.string_(self.secondary_mode)
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
@ -446,7 +446,7 @@ class ThermalScattering(EqualityMixin):
temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs]
table = cls(name, atomic_weight_ratio, kTs)
table.zaids = group.attrs['zaids']
table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']]
table.secondary_mode = group.attrs['secondary_mode'].decode()
# Read thermal elastic scattering
@ -628,8 +628,10 @@ class ThermalScattering(EqualityMixin):
table.elastic_mu_out[temperatures[0]].shape = \
(n_energy, n_mu)
# Get relevant ZAIDs
pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int)
table.zaids = pairs[np.nonzero(pairs)]
# Get relevant nuclides
for zaid, awr in ace.pairs:
if zaid > 0:
Z, A = divmod(zaid, 1000)
table.nuclides.append(ATOMIC_SYMBOL[Z] + str(A))
return table

21
openmc/mesh.py
View file

@ -187,15 +187,20 @@ class Mesh(object):
of the mesh.
For example the following code:
for mesh_index in mymesh.cell_generator():
print mesh_index
will produce the following output for a 3-D 2x2x2 mesh in mymesh:
[1, 1, 1]
[1, 1, 2]
[1, 2, 1]
[1, 2, 2]
...
.. code-block:: python
for mesh_index in mymesh.cell_generator():
print mesh_index
will produce the following output for a 3-D 2x2x2 mesh in mymesh::
[1, 1, 1]
[1, 1, 2]
[1, 2, 1]
[1, 2, 2]
...
"""

14
openmc/mgxs/mgxs.py
View file

@ -337,7 +337,7 @@ class MGXS(object):
if self.by_nuclide:
return self.get_nuclides()
else:
return 'sum'
return ['sum']
@property
def loaded_sp(self):
@ -1483,25 +1483,27 @@ class MGXS(object):
if self.by_nuclide and nuclides == 'sum':
# Use tally summation to sum across all nuclides
query_nuclides = self.get_nuclides()
xs_tally = self.xs_tally.summation(nuclides=query_nuclides)
query_nuclides = [nuclides]
xs_tally = self.xs_tally.summation(nuclides=self.get_nuclides())
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
# Remove nuclide column since it is homogeneous and redundant
if self.domain_type == 'mesh':
df.drop('nuclide', axis=1, level=0, inplace=True)
df.drop('sum(nuclide)', axis=1, level=0, inplace=True)
else:
df.drop('nuclide', axis=1, inplace=True)
df.drop('sum(nuclide)', axis=1, inplace=True)
# If the user requested a specific set of nuclides
elif self.by_nuclide and nuclides != 'all':
query_nuclides = nuclides
xs_tally = self.xs_tally.get_slice(nuclides=nuclides)
df = xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
# If the user requested all nuclides, keep nuclide column in dataframe
else:
query_nuclides = self.nuclides
df = self.xs_tally.get_pandas_dataframe(
distribcell_paths=distribcell_paths)
@ -1513,7 +1515,7 @@ class MGXS(object):
# Override energy groups bounds with indices
all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
all_groups = np.repeat(all_groups, self.num_nuclides)
all_groups = np.repeat(all_groups, len(query_nuclides))
if 'energy low [MeV]' in df and 'energyout low [MeV]' in df:
df.rename(columns={'energy low [MeV]': 'group in'},
inplace=True)

23
openmc/mgxs_library.py
View file

@ -99,9 +99,6 @@ class XSdata(object):
Unique identifier for the xsdata object
alias : str
Separate unique identifier for the xsdata object
zaid : int
1000*(atomic number) + mass number. As an example, the zaid of U235
would be 92235.
awr : float
Atomic weight ratio of an isotope. That is, the ratio of the mass
of the isotope to the mass of a single neutron.
@ -227,7 +224,6 @@ class XSdata(object):
self._energy_groups = energy_groups
self._representation = representation
self._alias = None
self._zaid = None
self._awr = None
self._kT = None
self._fissionable = False
@ -262,10 +258,6 @@ class XSdata(object):
def alias(self):
return self._alias
@property
def zaid(self):
return self._zaid
@property
def awr(self):
return self._awr
@ -396,13 +388,6 @@ class XSdata(object):
else:
self._alias = self._name
@zaid.setter
def zaid(self, zaid):
# Check type and value
check_type('zaid', zaid, Integral)
check_greater_than('zaid', zaid, 0)
self._zaid = zaid
@awr.setter
def awr(self, awr):
# Check validity of type and that the awr value is > 0
@ -1013,18 +998,10 @@ class XSdata(object):
subelement = ET.SubElement(element, 'kT')
subelement.text = str(self._kT)
if self._zaid is not None:
subelement = ET.SubElement(element, 'zaid')
subelement.text = str(self._zaid)
if self._awr is not None:
subelement = ET.SubElement(element, 'awr')
subelement.text = str(self._awr)
if self._kT is not None:
subelement = ET.SubElement(element, 'kT')
subelement.text = str(self._kT)
if self._fissionable is not None:
subelement = ET.SubElement(element, 'fissionable')
subelement.text = str(self._fissionable)

888
openmc/model/triso.py
View file

@ -1,13 +1,26 @@
from __future__ import division
import copy
from collections import Iterable
from numbers import Real
import warnings
import itertools
import random
from collections import Iterable, defaultdict
from numbers import Real
from random import uniform, gauss
from heapq import heappush, heappop
from math import pi, sin, cos, floor, log10, sqrt
from abc import ABCMeta, abstractproperty, abstractmethod
import numpy as np
try:
import scipy.spatial
_SCIPY_AVAILABLE = True
except ImportError:
_SCIPY_AVAILABLE = False
import openmc
import openmc.checkvalue as cv
class TRISO(openmc.Cell):
"""Tristructural-isotopic (TRISO) micro fuel particle
@ -82,6 +95,377 @@ class TRISO(openmc.Cell):
k_min:k_max+1, j_min:j_max+1, i_min:i_max+1]))
class _Domain(object):
"""Container in which to pack particles.
Parameters
----------
particle_radius : float
Radius of particles to be packed in container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
volume : float
Volume of the container.
"""
__metaclass__ = ABCMeta
def __init__(self, particle_radius, center=[0., 0., 0.]):
self._cell_length = None
self._limits = None
self.particle_radius = particle_radius
self.center = center
@property
def particle_radius(self):
return self._particle_radius
@property
def center(self):
return self._center
@abstractproperty
def limits(self):
pass
@abstractproperty
def cell_length(self):
pass
@abstractproperty
def volume(self):
pass
@particle_radius.setter
def particle_radius(self, particle_radius):
self._particle_radius = float(particle_radius)
self._limits = None
self._cell_length = None
@center.setter
def center(self, center):
if np.asarray(center).size != 3:
raise ValueError('Unable to set domain center to {} since it must '
'be of length 3'.format(center))
self._center = [float(x) for x in center]
self._limits = None
self._cell_length = None
def mesh_cell(self, p):
"""Calculate the index of the cell in a mesh overlaid on the domain in
which the given particle center falls.
Parameters
----------
p : Iterable of float
Cartesian coordinates of particle center.
Returns
-------
tuple of int
Indices of mesh cell.
"""
return tuple(int(p[i]/self.cell_length[i]) for i in range(3))
def nearby_mesh_cells(self, p):
"""Calculates the indices of all cells in a mesh overlaid on the domain
within one diameter of the given particle.
Parameters
----------
p : Iterable of float
Cartesian coordinates of particle center.
Returns
-------
list of tuple of int
Indices of mesh cells.
"""
d = 2*self.particle_radius
r = [[a/self.cell_length[i] for a in [p[i]-d, p[i], p[i]+d]]
for i in range(3)]
return list(itertools.product(*({int(x) for x in y} for y in r)))
@abstractmethod
def random_point(self):
"""Generate Cartesian coordinates of center of a particle that is
contained entirely within the domain with uniform probability.
Returns
-------
list of float
Cartesian coordinates of particle center.
"""
pass
class _CubicDomain(_Domain):
"""Cubic container in which to pack particles.
Parameters
----------
length : float
Length of each side of the cubic container.
particle_radius : float
Radius of particles to be packed in container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
length : float
Length of each side of the cubic container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, length, particle_radius, center=[0., 0., 0.]):
super(_CubicDomain, self).__init__(particle_radius, center)
self.length = length
@property
def length(self):
return self._length
@property
def limits(self):
if self._limits is None:
xlim = self.length/2 - self.particle_radius
self._limits = [[x - xlim for x in self.center],
[x + xlim for x in self.center]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [self.length, self.length, self.length]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return self.length**3
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
def random_point(self):
return [uniform(self.limits[0][0], self.limits[1][0]),
uniform(self.limits[0][1], self.limits[1][1]),
uniform(self.limits[0][2], self.limits[1][2])]
class _CylindricalDomain(_Domain):
"""Cylindrical container in which to pack particles.
Parameters
----------
length : float
Length along z-axis of the cylindrical container.
radius : float
Radius of the cylindrical container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
length : float
Length along z-axis of the cylindrical container.
radius : float
Radius of the cylindrical container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]):
super(_CylindricalDomain, self).__init__(particle_radius, center)
self.length = length
self.radius = radius
@property
def length(self):
return self._length
@property
def radius(self):
return self._radius
@property
def limits(self):
if self._limits is None:
xlim = self.length/2 - self.particle_radius
rlim = self.radius - self.particle_radius
self._limits = [[self.center[0] - rlim, self.center[1] - rlim,
self.center[2] - xlim],
[self.center[0] + rlim, self.center[1] + rlim,
self.center[2] + xlim]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [2*self.radius, 2*self.radius, self.length]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return self.length * pi * self.radius**2
@length.setter
def length(self, length):
self._length = float(length)
self._limits = None
self._cell_length = None
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
def random_point(self):
r = sqrt(uniform(0, (self.radius - self.particle_radius)**2))
t = uniform(0, 2*pi)
return [r*cos(t) + self.center[0], r*sin(t) + self.center[1],
uniform(self.limits[0][2], self.limits[1][2])]
class _SphericalDomain(_Domain):
"""Spherical container in which to pack particles.
Parameters
----------
radius : float
Radius of the spherical container.
center : Iterable of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
Attributes
----------
radius : float
Radius of the spherical container.
particle_radius : float
Radius of particles to be packed in container.
center : list of float
Cartesian coordinates of the center of the container. Default is
[0., 0., 0.]
cell_length : list of float
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
volume : float
Volume of the container.
"""
def __init__(self, radius, particle_radius, center=[0., 0., 0.]):
super(_SphericalDomain, self).__init__(particle_radius, center)
self.radius = radius
@property
def radius(self):
return self._radius
@property
def limits(self):
if self._limits is None:
rlim = self.radius - self.particle_radius
self._limits = [[x - rlim for x in self.center],
[x + rlim for x in self.center]]
return self._limits
@property
def cell_length(self):
if self._cell_length is None:
mesh_length = [2*self.radius, 2*self.radius, 2*self.radius]
self._cell_length = [x/int(x/(4*self.particle_radius))
for x in mesh_length]
return self._cell_length
@property
def volume(self):
return 4/3 * pi * self.radius**3
@radius.setter
def radius(self, radius):
self._radius = float(radius)
self._limits = None
self._cell_length = None
@limits.setter
def limits(self, limits):
self._limits = limits
def random_point(self):
x = (gauss(0, 1), gauss(0, 1), gauss(0, 1))
r = (uniform(0, (self.radius - self.particle_radius)**3)**(1/3) /
sqrt(x[0]**2 + x[1]**2 + x[2]**2))
return [r*x[i] + self.center[i] for i in range(3)]
def create_triso_lattice(trisos, lower_left, pitch, shape, background):
"""Create a lattice containing TRISO particles for optimized tracking.
@ -153,3 +537,503 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background):
lattice.outer = openmc.Universe(cells=[background_cell])
return lattice
def _random_sequential_pack(domain, n_particles):
"""Random sequential packing of particles within a container.
Parameters
----------
domain : openmc.model._Domain
Container in which to pack particles.
n_particles : int
Number of particles to pack.
Returns
------
numpy.ndarray
Cartesian coordinates of centers of particles.
"""
sqd = (2*domain.particle_radius)**2
particles = []
mesh = defaultdict(list)
for i in range(n_particles):
# Randomly sample new center coordinates while there are any overlaps
while True:
p = domain.random_point()
idx = domain.mesh_cell(p)
if any((p[0]-q[0])**2 + (p[1]-q[1])**2 + (p[2]-q[2])**2 < sqd
for q in mesh[idx]):
continue
else:
break
particles.append(p)
for idx in domain.nearby_mesh_cells(p):
mesh[idx].append(p)
return np.array(particles)
def _close_random_pack(domain, particles, contraction_rate):
"""Close random packing of particles using the Jodrey-Tory algorithm.
Parameters
----------
domain : openmc.model._Domain
Container in which to pack particles.
particles : numpy.ndarray
Initial Cartesian coordinates of centers of particles.
contraction_rate : float
Contraction rate of outer diameter.
"""
def add_rod(d, i, j):
"""Add a new rod to the priority queue.
Parameters
----------
d : float
distance between centers of particles i and j.
i, j : int
Index of particles in particles array.
"""
rod = [d, i, j]
rods_map[i] = (j, rod)
rods_map[j] = (i, rod)
heappush(rods, rod)
def remove_rod(i):
"""Mark the rod containing particle i as removed.
Parameters
----------
i : int
Index of particle in particles array.
"""
if i in rods_map:
j, rod = rods_map.pop(i)
del rods_map[j]
rod[1] = removed
rod[2] = removed
def pop_rod():
"""Remove and return the shortest rod.
Returns
-------
d : float
distance between centers of particles i and j.
i, j : int
Index of particles in particles array.
"""
while rods:
d, i, j = heappop(rods)
if i != removed and j != removed:
del rods_map[i]
del rods_map[j]
return d, i, j
def create_rod_list():
"""Generate sorted list of rods (distances between particle centers).
Rods are arranged in a heap where each element contains the rod length
and the particle indices. A rod between particles p and q is only
included if the distance between p and q could not be changed by the
elimination of a greater overlap, i.e. q has no nearer neighbors than p.
A mapping of particle ids to rods is maintained in 'rods_map'. Each key
in the dict is the id of a particle that is in the rod list, and the
value is the id of its nearest neighbor and the rod that contains them.
The dict is used to find rods in the priority queue and to mark removed
rods so rods can be "removed" without breaking the heap structure
invariant.
"""
# Create KD tree for quick nearest neighbor search
tree = scipy.spatial.cKDTree(particles)
# Find distance to nearest neighbor and index of nearest neighbor for
# all particles
d, n = tree.query(particles, k=2)
d = d[:,1]
n = n[:,1]
# Array of particle indices, indices of nearest neighbors, and
# distances to nearest neighbors
a = np.vstack((list(range(n.size)), n, d)).T
# Sort along second column and swap first and second columns to create
# array of nearest neighbor indices, indices of particles they are
# nearest neighbors of, and distances between them
b = a[a[:,1].argsort()]
b[:,[0, 1]] = b[:,[1, 0]]
# Find the intersection between 'a' and 'b': a list of particles who
# are each other's nearest neighbors and the distance between them
r = list({tuple(x) for x in a} & {tuple(x) for x in b})
# Remove duplicate rods and sort by distance
r = map(list, set([(x[2], int(min(x[0:2])), int(max(x[0:2])))
for x in r]))
# Clear priority queue and add rods
del rods[:]
rods_map.clear()
for d, i, j in r:
add_rod(d, i, j)
# Inner diameter is set initially to the shortest center-to-center
# distance between any two particles
if rods:
inner_diameter[0] = rods[0][0]
def update_mesh(i):
"""Update which mesh cells the particle is in based on new particle
center coordinates.
'mesh'/'mesh_map' is a two way dictionary used to look up which
particles are located within one diameter of a given mesh cell and
which mesh cells a given particle center is within one diameter of.
This is used to speed up the nearest neighbor search.
Parameters
----------
i : int
Index of particle in particles array.
"""
# Determine which mesh cells the particle is in and remove the
# particle id from those cells
for idx in mesh_map[i]:
mesh[idx].remove(i)
del mesh_map[i]
# Determine which mesh cells are within one diameter of particle's
# center and add this particle to the list of particles in those cells
for idx in domain.nearby_mesh_cells(particles[i]):
mesh[idx].add(i)
mesh_map[i].add(idx)
def reduce_outer_diameter():
"""Reduce the outer diameter so that at the (i+1)-st iteration it is:
d_out^(i+1) = d_out^(i) - (1/2)^(j) * d_out0 * k / n,
where k is the contraction rate, n is the number of particles, and
j = floor(-log10(pf_out - pf_in)).
"""
inner_pf = (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles /
domain.volume)
outer_pf = (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles /
domain.volume)
j = floor(-log10(outer_pf - inner_pf))
outer_diameter[0] = (outer_diameter[0] - 0.5**j * contraction_rate *
initial_outer_diameter / n_particles)
def repel_particles(i, j, d):
"""Move particles p and q apart according to the following
transformation (accounting for reflective boundary conditions on
domain):
r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n))
r_j^(n+1) = r_j^(n) - 1/2(d_out^(n+1) - d^(n))
Parameters
----------
i, j : int
Index of particles in particles array.
d : float
distance between centers of particles i and j.
"""
# Moving each particle distance 'r' away from the other along the line
# joining the particle centers will ensure their final distance is equal
# to the outer diameter
r = (outer_diameter[0] - d)/2
v = (particles[i] - particles[j])/d
particles[i] += r*v
particles[j] -= r*v
# Apply reflective boundary conditions
particles[i] = particles[i].clip(domain.limits[0], domain.limits[1])
particles[j] = particles[j].clip(domain.limits[0], domain.limits[1])
update_mesh(i)
update_mesh(j)
def nearest(i):
"""Find index of nearest neighbor of particle i.
Parameters
----------
i : int
Index in particles array of particle for which to find nearest
neighbor.
Returns
-------
int
Index in particles array of nearest neighbor of i
float
distance between i and nearest neighbor.
"""
# Need the second nearest neighbor of i since the nearest neighbor
# will be itself. Using argpartition, the k-th nearest neighbor is
# placed at index k.
idx = list(mesh[domain.mesh_cell(particles[i])])
dists = scipy.spatial.distance.cdist([particles[i]], particles[idx])[0]
if dists.size > 1:
j = dists.argpartition(1)[1]
return idx[j], dists[j]
else:
return None, None
def update_rod_list(i, j):
"""Update the rod list with the new nearest neighbors of particles i
and j since their overlap was eliminated.
Parameters
----------
i, j : int
Index of particles in particles array.
"""
# If the nearest neighbor k of particle i has no nearer neighbors,
# remove the rod currently containing k from the rod list and add rod
# k-i, keeping the rod list sorted
k, d_ik = nearest(i)
if k and nearest(k)[0] == i:
remove_rod(k)
add_rod(d_ik, i, k)
l, d_jl = nearest(j)
if l and nearest(l)[0] == j:
remove_rod(l)
add_rod(d_jl, j, l)
# Set inner diameter to the shortest distance between two particle
# centers
if rods:
inner_diameter[0] = rods[0][0]
if not _SCIPY_AVAILABLE:
raise ImportError('SciPy must be installed to perform '
'close random packing.')
n_particles = len(particles)
diameter = 2*domain.particle_radius
# Flag for marking rods that have been removed from priority queue
removed = -1
# Outer diameter initially set to arbitrary value that yields pf of 1
initial_outer_diameter = 2*(domain.volume/(n_particles*4/3*pi))**(1/3)
# Inner and outer diameter of particles will change during packing
outer_diameter = [initial_outer_diameter]
inner_diameter = [0]
rods = []
rods_map = {}
mesh = defaultdict(set)
mesh_map = defaultdict(set)
for i in range(n_particles):
for idx in domain.nearby_mesh_cells(particles[i]):
mesh[idx].add(i)
mesh_map[i].add(idx)
while True:
create_rod_list()
if inner_diameter[0] >= diameter:
break
while True:
d, i, j = pop_rod()
reduce_outer_diameter()
repel_particles(i, j, d)
update_rod_list(i, j)
if inner_diameter[0] >= diameter or not rods:
break
def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
domain_radius=None, domain_center=[0., 0., 0.],
n_particles=None, packing_fraction=None,
initial_packing_fraction=0.3, contraction_rate=1/400, seed=1):
"""Generate a random, non-overlapping configuration of TRISO particles
within a container.
Parameters
----------
radius : float
Outer radius of TRISO particles.
fill : openmc.Universe
Universe which contains all layers of the TRISO particle.
domain_shape : {'cube', 'cylinder', or 'sphere'}
Geometry of the container in which the TRISO particles are packed.
domain_length : float
Length of the container (if cube or cylinder).
domain_radius : float
Radius of the container (if cylinder or sphere).
domain_center : Iterable of float
Cartesian coordinates of the center of the container.
n_particles : int
Number of TRISO particles to pack in the domain. Exactly one of
'n_particles' and 'packing_fraction' should be specified -- the other
will be calculated.
packing_fraction : float
Packing fraction of particles. Exactly one of 'n_particles' and
'packing_fraction' should be specified -- the other will be calculated.
initial_packing_fraction : float, optional
Packing fraction used to initialize the configuration of particles in
the domain. Default value is 0.3. It is not recommended to set the
initial packing fraction much higher than 0.3 as the random sequential
packing algorithm becomes prohibitively slow as it approaches its limit
(~0.38).
contraction_rate : float, optional
Contraction rate of outer diameter. This can affect the speed of the
close random packing algorithm. Default value is 1/400.
seed : int, optional
RNG seed.
Returns
-------
trisos : list of openmc.model.TRISO
List of TRISO particles in the domain.
Notes
-----
The particle configuration is generated using a combination of random
sequential packing (RSP) and close random packing (CRP). RSP performs
better than CRP for lower packing fractions (pf), but it becomes
prohibitively slow as it approaches its packing limit (~0.38). CRP can
achieve higher pf of up to ~0.64 and scales better with increasing pf.
If the desired pf is below some threshold for which RSP will be faster than
CRP ('initial_packing_fraction'), only RSP is used. If a higher pf is
required, particles with a radius smaller than the desired final radius
(and therefore with a smaller pf) are initialized within the domain using
RSP. This initial configuration of particles is then used as a starting
point for CRP using Jodrey and Tory's algorithm [1]_.
In RSP, particle centers are placed one by one at random, and placement
attempts for a particle are made until the particle is not overlapping any
others. This implementation of the algorithm uses a mesh over the domain
to speed up the nearest neighbor search by only searching for a particle's
neighbors within that mesh cell.
In CRP, each particle is assigned two diameters, and inner and an outer,
which approach each other during the simulation. The inner diameter,
defined as the minimum center-to-center distance, is the true diameter of
the particles and defines the pf. At each iteration the worst overlap
between particles based on outer diameter is eliminated by moving the
particles apart along the line joining their centers. Iterations continue
until the two diameters converge or until the desired pf is reached.
References
----------
.. [1] W. S. Jodrey and E. M. Tory, "Computer simulation of close random
packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351.
"""
# Check for valid container geometry and dimensions
if domain_shape not in ['cube', 'cylinder', 'sphere']:
raise ValueError('Unable to set domain_shape to "{}". Only "cube", '
'"cylinder", and "sphere" are '
'supported."'.format(domain_shape))
if not domain_length and domain_shape in ['cube', 'cylinder']:
raise ValueError('"domain_length" must be specified for {} domain '
'geometry '.format(domain_shape))
if not domain_radius and domain_shape in ['cylinder', 'sphere']:
raise ValueError('"domain_radius" must be specified for {} domain '
'geometry '.format(domain_shape))
if domain_shape is 'cube':
domain = _CubicDomain(length=domain_length, particle_radius=radius,
center=domain_center)
elif domain_shape is 'cylinder':
domain = _CylindricalDomain(length=domain_length, radius=domain_radius,
particle_radius=radius, center=domain_center)
elif domain_shape is 'sphere':
domain = _SphericalDomain(radius=domain_radius, particle_radius=radius,
center=domain_center)
# Calculate the packing fraction if the number of particles is specified;
# otherwise, calculate the number of particles from the packing fraction.
if ((n_particles is None and packing_fraction is None) or
(n_particles is not None and packing_fraction is not None)):
raise ValueError('Exactly one of "n_particles" and "packing_fraction" '
'must be specified.')
elif packing_fraction is None:
n_particles = int(n_particles)
packing_fraction = 4/3*pi*radius**3*n_particles / domain.volume
elif n_particles is None:
packing_fraction = float(packing_fraction)
n_particles = int(packing_fraction*domain.volume // (4/3*pi*radius**3))
# Check for valid packing fractions for each algorithm
if packing_fraction >= 0.64:
raise ValueError('Packing fraction of {} is greater than the '
'packing fraction limit for close random '
'packing (0.64)'.format(packing_fraction))
if initial_packing_fraction >= 0.38:
raise ValueError('Initial packing fraction of {} is greater than the '
'packing fraction limit for random sequential'
'packing (0.38)'.format(initial_packing_fraction))
if initial_packing_fraction > packing_fraction:
initial_packing_fraction = packing_fraction
if packing_fraction > 0.3:
initial_packing_fraction = 0.3
random.seed(seed)
# Calculate the particle radius used in the initial random sequential
# packing from the initial packing fraction
initial_radius = (3/4 * initial_packing_fraction * domain.volume /
(pi * n_particles))**(1/3)
domain.particle_radius = initial_radius
# Recalculate the limits for the initial random sequential packing using
# the desired final particle radius to ensure particles are fully contained
# within the domain during the close random pack
domain.limits = [[x - initial_radius + radius for x in domain.limits[0]],
[x + initial_radius - radius for x in domain.limits[1]]]
# Generate non-overlapping particles for an initial inner radius using
# random sequential packing algorithm
particles = _random_sequential_pack(domain, n_particles)
# Use the particle configuration produced in random sequential packing as a
# starting point for close random pack with the desired final particle
# radius
if initial_packing_fraction != packing_fraction:
domain.particle_radius = radius
_close_random_pack(domain, particles, contraction_rate)
trisos = []
for p in particles:
trisos.append(TRISO(radius, fill, p))
return trisos

16
openmc/nuclide.py
View file

@ -20,9 +20,6 @@ class Nuclide(object):
----------
name : str
Name of the nuclide, e.g. U235
zaid : int
1000*(atomic number) + mass number. As an example, the zaid of U235
would be 92235.
scattering : 'data' or 'iso-in-lab' or None
The type of angular scattering distribution to use
@ -31,7 +28,6 @@ class Nuclide(object):
def __init__(self, name=''):
# Initialize class attributes
self._name = ''
self._zaid = None
self._scattering = None
# Set the Material class attributes
@ -62,8 +58,6 @@ class Nuclide(object):
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)
if self.zaid is not None:
string += '{0: <16}{1}{2}\n'.format('\tZAID', '=\t', self.zaid)
if self.scattering is not None:
string += '{0: <16}{1}{2}\n'.format('\tscattering', '=\t',
self.scattering)
@ -73,10 +67,6 @@ class Nuclide(object):
def name(self):
return self._name
@property
def zaid(self):
return self._zaid
@property
def scattering(self):
return self._scattering
@ -96,14 +86,8 @@ class Nuclide(object):
'"{}" is being renamed as "{}".'.format(name, self._name)
warnings.warn(msg)
@zaid.setter
def zaid(self, zaid):
check_type('zaid', zaid, Integral)
self._zaid = zaid
@scattering.setter
def scattering(self, scattering):
if not scattering in ['data', 'iso-in-lab']:
msg = 'Unable to set scattering for Nuclide to {0} ' \
'which is not "data" or "iso-in-lab"'.format(scattering)

86
openmc/plots.py
View file

@ -67,6 +67,11 @@ class Plot(object):
col_spec : dict
Dictionary indicating that certain cells/materials (keys) should be
colored with a specific RGB (values)
level : int
Universe depth to plot at
meshlines : dict
Dictionary defining type, id, linewidth and color of a regular mesh
to be plotted on top of a plot
"""
@ -81,10 +86,12 @@ class Plot(object):
self._color = 'cell'
self._type = 'slice'
self._basis = 'xy'
self._background = [0, 0, 0]
self._background = None
self._mask_components = None
self._mask_background = None
self._col_spec = None
self._level = None
self._meshlines = None
@property
def id(self):
@ -138,6 +145,14 @@ class Plot(object):
def col_spec(self):
return self._col_spec
@property
def level(self):
return self._level
@property
def meshlines(self):
return self._meshlines
@id.setter
def id(self, plot_id):
if plot_id is None:
@ -231,9 +246,9 @@ class Plot(object):
@mask_components.setter
def mask_components(self, mask_components):
cv.check_type('plot mask_components', mask_components, Iterable, Integral)
cv.check_type('plot mask components', mask_components, Iterable, Integral)
for component in mask_components:
cv.check_greater_than('plot mask_components', component, 0, True)
cv.check_greater_than('plot mask components', component, 0, True)
self._mask_components = mask_components
@mask_background.setter
@ -245,6 +260,45 @@ class Plot(object):
cv.check_less_than('plot mask background', rgb, 256)
self._mask_background = mask_background
@level.setter
def level(self, plot_level):
cv.check_type('plot level', plot_level, Integral)
cv.check_greater_than('plot level', plot_level, 0, equality=True)
self._level = plot_level
@meshlines.setter
def meshlines(self, meshlines):
cv.check_type('plot meshlines', meshlines, dict)
if 'type' not in meshlines:
msg = 'Unable to set on plot the meshlines "{0}" which ' \
'does not have a "type" key'.format(meshlines)
raise ValueError(msg)
elif meshlines['type'] not in ['tally', 'entropy', 'ufs', 'cmfd']:
msg = 'Unable to set the meshlines with ' \
'type "{0}"'.format(meshlines['type'])
raise ValueError(msg)
if 'id' in meshlines:
cv.check_type('plot meshlines id', meshlines['id'], Integral)
cv.check_greater_than('plot meshlines id', meshlines['id'], 0,
equality=True)
if 'linewidth' in meshlines:
cv.check_type('plot mesh linewidth', meshlines['linewidth'], Integral)
cv.check_greater_than('plot mesh linewidth', meshlines['linewidth'],
0, equality=True)
if 'color' in meshlines:
cv.check_type('plot meshlines color', meshlines['color'], Iterable,
Integral)
cv.check_length('plot meshlines color', meshlines['color'], 3)
for rgb in meshlines['color']:
cv.check_greater_than('plot meshlines color', rgb, 0, True)
cv.check_less_than('plot meshlines color', rgb, 256)
self._meshlines = meshlines
def __repr__(self):
string = 'Plot\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -256,11 +310,16 @@ class Plot(object):
string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin)
string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
string += '{0: <16}{1}{2}\n'.format('\tBackground', '=\t',
self._background)
string += '{0: <16}{1}{2}\n'.format('\tMask components', '=\t',
self._mask_components)
string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t',
string += '{0: <16}{1}{2}\n'.format('\tMask background', '=\t',
self._mask_background)
string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
string += '{0: <16}{1}{2}\n'.format('\tLevel', '=\t', self._level)
string += '{0: <16}{1}{2}\n'.format('\tMeshlines', '=\t',
self._meshlines)
return string
def colorize(self, geometry, seed=1):
@ -382,7 +441,7 @@ class Plot(object):
subelement = ET.SubElement(element, "pixels")
subelement.text = ' '.join(map(str, self._pixels))
if self._mask_background is not None:
if self._background is not None:
subelement = ET.SubElement(element, "background")
subelement.text = ' '.join(map(str, self._background))
@ -400,6 +459,21 @@ class Plot(object):
subelement.set("background", ' '.join(map(
str, self._mask_background)))
if self._level is not None:
subelement = ET.SubElement(element, "level")
subelement.text = str(self._level)
if self._meshlines is not None:
subelement = ET.SubElement(element, "meshlines")
subelement.set("meshtype", self._meshlines['type'])
if self._meshlines['id'] is not None:
subelement.set("id", str(self._meshlines['id']))
if self._meshlines['linewidth'] is not None:
subelement.set("linewidth", str(self._meshlines['linewidth']))
if self._meshlines['color'] is not None:
subelement.set("color", ' '.join(map(
str, self._meshlines['color'])))
return element

122
openmc/settings.py
View file

@ -1,4 +1,4 @@
from collections import Iterable, MutableSequence
from collections import Iterable, MutableSequence, Mapping
from numbers import Real, Integral
import warnings
from xml.etree import ElementTree as ET
@ -78,8 +78,6 @@ class Settings(object):
cross section library. If it is not set, the
:envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. A
multipole library is optional.
energy_grid : {'nuclide', 'logarithm', 'material-union'}
Set the method used to search energy grids.
energy_mode : {'continuous-energy', 'multi-group'}
Set whether the calculation should be continuous-energy or multi-group.
max_order : int
@ -103,6 +101,14 @@ class Settings(object):
Coordinates of the lower-left point of the Shannon entropy mesh
entropy_upper_right : tuple or list
Coordinates of the upper-right point of the Shannon entropy mesh
temperature : dict
Defines a default temperature and method for treating intermediate
temperatures at which nuclear data doesn't exist. Accepted keys are
'default', 'method', and 'tolerance'. The value for 'default' should be
a float representing the default temperature in Kelvin. The value for
'method' should be 'nearest' or 'interpolation'. If the method is
'nearest', 'tolerance' indicates a range of temperature within which
cross sections may be used.
trigger_active : bool
Indicate whether tally triggers are used
trigger_max_batches : int
@ -130,9 +136,6 @@ class Settings(object):
Coordinates of the lower-left point of the UFS mesh
ufs_upper_right : tuple or list
Coordinates of the upper-right point of the UFS mesh
use_windowed_multipole : bool
Whether or not windowed multipole can be used to evaluate resolved
resonance cross sections.
resonance_scattering : ResonanceScattering or iterable of ResonanceScattering
The elastic scattering model to use for resonant isotopes
volume_calculations : VolumeCalculation or iterable of VolumeCalculation
@ -160,7 +163,6 @@ class Settings(object):
self._confidence_intervals = None
self._cross_sections = None
self._multipole_library = None
self._energy_grid = None
self._ptables = None
self._run_cmfd = None
self._seed = None
@ -197,6 +199,8 @@ class Settings(object):
self._trace = None
self._track = None
self._temperature = {}
# Cutoff subelement
self._weight = None
self._weight_avg = None
@ -216,7 +220,6 @@ class Settings(object):
self._settings_file = ET.Element("settings")
self._run_mode_subelement = None
self._multipole_active = None
self._resonance_scattering = cv.CheckedList(
ResonanceScattering, 'resonance scattering models')
@ -271,10 +274,6 @@ class Settings(object):
def multipole_library(self):
return self._multipole_library
@property
def energy_grid(self):
return self._energy_grid
@property
def ptables(self):
return self._ptables
@ -363,6 +362,10 @@ class Settings(object):
def verbosity(self):
return self._verbosity
@property
def temperature(self):
return self._temperature
@property
def trace(self):
return self._trace
@ -415,10 +418,6 @@ class Settings(object):
def dd_count_interactions(self):
return self._dd_count_interactions
@property
def use_windowed_multipole(self):
return self._multipole_active
@property
def resonance_scattering(self):
return self._resonance_scattering
@ -593,12 +592,6 @@ class Settings(object):
cv.check_type('cross sections', multipole_library, basestring)
self._multipole_library = multipole_library
@energy_grid.setter
def energy_grid(self, energy_grid):
cv.check_value('energy grid', energy_grid,
['nuclide', 'logarithm', 'material-union'])
self._energy_grid = energy_grid
@ptables.setter
def ptables(self, ptables):
cv.check_type('probability tables', ptables, bool)
@ -674,6 +667,21 @@ class Settings(object):
cv.check_type('no reduction option', no_reduce, bool)
self._no_reduce = no_reduce
@temperature.setter
def temperature(self, temperature):
cv.check_type('temperature settings', temperature, Mapping)
for key, value in temperature.items():
cv.check_value('temperature key', key,
['default', 'method', 'tolerance'])
if key == 'default':
cv.check_type('default temperature', value, Real)
elif key == 'method':
cv.check_value('temperature method', value,
['nearest', 'interpolation', 'multipole'])
elif key == 'tolerance':
cv.check_type('temperature tolerance', value, Real)
self._temperature = temperature
@threads.setter
def threads(self, threads):
cv.check_type('number of threads', threads, Integral)
@ -801,11 +809,6 @@ class Settings(object):
self._dd_count_interactions = interactions
@use_windowed_multipole.setter
def use_windowed_multipole(self, active):
cv.check_type('use_windowed_multipole', active, bool)
self._multipole_active = active
@resonance_scattering.setter
def resonance_scattering(self, res):
if not isinstance(res, MutableSequence):
@ -963,11 +966,6 @@ class Settings(object):
element = ET.SubElement(self._settings_file, "multipole_library")
element.text = str(self._multipole_library)
def _create_energy_grid_subelement(self):
if self._energy_grid is not None:
element = ET.SubElement(self._settings_file, "energy_grid")
element.text = str(self._energy_grid)
def _create_ptables_subelement(self):
if self._ptables is not None:
element = ET.SubElement(self._settings_file, "ptables")
@ -1050,6 +1048,13 @@ class Settings(object):
element = ET.SubElement(self._settings_file, "no_reduce")
element.text = str(self._no_reduce).lower()
def _create_temperature_subelements(self):
if self.temperature:
for key, value in self.temperature.items():
element = ET.SubElement(self._settings_file,
"temperature_{}".format(key))
element.text = str(value)
def _create_threads_subelement(self):
if self._threads is not None:
element = ET.SubElement(self._settings_file, "threads")
@ -1107,20 +1112,10 @@ class Settings(object):
subelement = ET.SubElement(element, "count_interactions")
subelement.text = str(self._dd_count_interactions).lower()
def _create_use_multipole_subelement(self):
if self._multipole_active is not None:
element = ET.SubElement(self._settings_file,
"use_windowed_multipole")
element.text = str(self._multipole_active)
def _create_resonance_scattering_subelement(self):
if len(self.resonance_scattering) > 0:
elem = ET.SubElement(self._settings_file, 'resonance_scattering')
for r in self.resonance_scattering:
if r.nuclide.name != r.nuclide_0K.name:
raise ValueError("The nuclide and nuclide_0K attributes of "
"a ResonantScattering object must have "
"identical names.")
elem.append(r.to_xml_element())
def export_to_xml(self):
@ -1142,7 +1137,6 @@ class Settings(object):
self._create_confidence_intervals()
self._create_cross_sections_subelement()
self._create_multipole_library_subelement()
self._create_energy_grid_subelement()
self._create_energy_mode_subelement()
self._create_max_order_subelement()
self._create_ptables_subelement()
@ -1155,11 +1149,11 @@ class Settings(object):
self._create_no_reduce_subelement()
self._create_threads_subelement()
self._create_verbosity_subelement()
self._create_temperature_subelements()
self._create_trace_subelement()
self._create_track_subelement()
self._create_ufs_subelement()
self._create_dd_subelement()
self._create_use_multipole_subelement()
self._create_resonance_scattering_subelement()
self._create_volume_calcs_subelement()
@ -1175,14 +1169,26 @@ class Settings(object):
class ResonanceScattering(object):
"""Specification of the elastic scattering model for resonant isotopes
Parameters
----------
nuclide : openmc.Nuclide
The nuclide affected by this resonance scattering treatment.
method : {'ARES', 'CXS', 'DBRC', 'WCM'}
The method used to sample outgoing scattering energies. Valid options
are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening
rejection correction), and 'WCM' (weight correction method).
E_min : float
The minimum energy above which the specified method is applied. By
default, CXS will be used below E_min.
E_max : float
The maximum energy below which the specified method is applied. By
default, the asymptotic target-at-rest model is applied above E_max.
Attributes
----------
nuclide : openmc.Nuclide
The nuclide affected by this resonance scattering treatment.
nuclide_0K : openmc.Nuclide
This should be the same isotope as the nuclide attribute above, but it
should have an xs attribute that identifies 0 Kelvin data.
method : str
method : {'ARES', 'CXS', 'DBRC', 'WCM'}
The method used to sample outgoing scattering energies. Valid options
are 'ARES', 'CXS' (constant cross section), 'DBRC' (Doppler broadening
rejection correction), and 'WCM' (weight correction method).
@ -1195,21 +1201,20 @@ class ResonanceScattering(object):
"""
def __init__(self):
self._nuclide = None
self._nuclide_0K = None
self._method = None
def __init__(self, nuclide, method='CXS', E_min=None, E_max=None):
self._E_min = None
self._E_max = None
self.nuclide = nuclide
self.method = method
if E_min is not None:
self.E_min = E_min
if E_max is not None:
self.E_max = E_max
@property
def nuclide(self):
return self._nuclide
@property
def nuclide_0K(self):
return self._nuclide_0K
@property
def method(self):
return self._method
@ -1227,11 +1232,6 @@ class ResonanceScattering(object):
cv.check_type('nuclide', nuc, Nuclide)
self._nuclide = nuc
@nuclide_0K.setter
def nuclide_0K(self, nuc):
cv.check_type('nuclide_0K', nuc, Nuclide)
self._nuclide_0K = nuc
@method.setter
def method(self, m):
cv.check_value('method', m, ('ARES', 'CXS', 'DBRC', 'WCM'))

4
openmc/summary.py
View file

@ -83,11 +83,10 @@ class Summary(object):
n_nuclides = self._f['nuclides/n_nuclides_total'].value
names = self._f['nuclides/names'].value
awrs = self._f['nuclides/awrs'].value
zaids = self._f['nuclides/zaids'].value
for n in range(n_nuclides):
name = names[n].decode()
name = name[:name.find('.')]
self.nuclides[name] = (zaids[n], awrs[n])
self.nuclides[name] = awrs[n]
def _read_geometry(self):
# Read in and initialize the Materials and Geometry
@ -113,7 +112,6 @@ class Summary(object):
material_id = int(key.lstrip('material '))
index = self._f['materials'][key]['index'].value
name = self._f['materials'][key]['name'].value.decode()
temperature = self._f['materials'][key]['temperature'].value.decode()
density = self._f['materials'][key]['atom_density'].value
nuc_densities = self._f['materials'][key]['nuclide_densities'][...]
nuclides = self._f['materials'][key]['nuclides'].value

275
src/algorithm.F90 Normal file
View file

@ -0,0 +1,275 @@
module algorithm
use constants
use stl_vector, only: VectorInt, VectorReal
implicit none
integer, parameter :: MAX_ITERATION = 64
interface binary_search
module procedure binary_search_real, binary_search_int4, binary_search_int8
end interface binary_search
interface sort
module procedure sort_int, sort_real, sort_vector_int, sort_vector_real
end interface sort
interface find
module procedure find_int, find_real, find_vector_int, find_vector_real
end interface find
contains
!===============================================================================
! BINARY_SEARCH performs a binary search of an array to find where a specific
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
pure function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
real(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_real
pure function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
integer, intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int4
pure function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
integer(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int8
!===============================================================================
! SORT sorts an array in place using an insertion sort.
!===============================================================================
pure subroutine sort_int(array)
integer, intent(inout) :: array(:)
integer :: k, m
integer :: temp
if (size(array) > 1) then
SORT: do k = 2, size(array)
! Save value to move
m = k
temp = array(k)
MOVE_OVER: do while (m > 1)
! Check if insertion value is greater than (m-1)th value
if (temp >= array(m - 1)) exit
! Move values over until hitting one that's not larger
array(m) = array(m - 1)
m = m - 1
end do MOVE_OVER
! Put the original value into its new position
array(m) = temp
end do SORT
end if
end subroutine sort_int
pure subroutine sort_real(array)
real(8), intent(inout) :: array(:)
integer :: k, m
real(8) :: temp
if (size(array) > 1) then
SORT: do k = 2, size(array)
! Save value to move
m = k
temp = array(k)
MOVE_OVER: do while (m > 1)
! Check if insertion value is greater than (m-1)th value
if (temp >= array(m - 1)) exit
! Move values over until hitting one that's not larger
array(m) = array(m - 1)
m = m - 1
end do MOVE_OVER
! Put the original value into its new position
array(m) = temp
end do SORT
end if
end subroutine sort_real
pure subroutine sort_vector_int(vec)
type(VectorInt), intent(inout) :: vec
call sort_int(vec % data(1:vec%size()))
end subroutine sort_vector_int
pure subroutine sort_vector_real(vec)
type(VectorReal), intent(inout) :: vec
call sort_real(vec % data(1:vec%size()))
end subroutine sort_vector_real
!===============================================================================
! FIND determines the index of the first occurrence of a value in an array. If
! the value does not appear in the array, -1 is returned.
!===============================================================================
pure function find_int(array, val) result(index)
integer, intent(in) :: array(:)
integer, intent(in) :: val
integer :: index
integer :: i
index = -1
do i = 1, size(array)
if (array(i) == val) then
index = i
exit
end if
end do
end function find_int
pure function find_real(array, val) result(index)
real(8), intent(in) :: array(:)
real(8), intent(in) :: val
integer :: index
integer :: i
index = -1
do i = 1, size(array)
if (array(i) == val) then
index = i
exit
end if
end do
end function find_real
pure function find_vector_int(vec, val) result(index)
type(VectorInt), intent(in) :: vec
integer, intent(in) :: val
integer :: index
index = find_int(vec % data(1:vec % size()), val)
end function find_vector_int
pure function find_vector_real(vec, val) result(index)
type(VectorReal), intent(in) :: vec
real(8), intent(in) :: val
integer :: index
index = find_real(vec % data(1:vec % size()), val)
end function find_vector_real
end module algorithm

2
src/angle_distribution.F90
View file

@ -2,12 +2,12 @@ module angle_distribution
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: read_attribute, get_shape, read_dataset, &
open_dataset, close_dataset
use random_lcg, only: prn
use search, only: binary_search
implicit none
private

2
src/cmfd_execute.F90
View file

@ -213,13 +213,13 @@ contains
subroutine cmfd_reweight(new_weights)
use algorithm, only: binary_search
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use global, only: meshes, source_bank, work, n_user_meshes, cmfd, &
master
use mesh_header, only: RegularMesh
use mesh, only: count_bank_sites, get_mesh_indices
use search, only: binary_search
use string, only: to_str
#ifdef MPI

12
src/constants.F90
View file

@ -268,6 +268,12 @@ module constants
JENDL_33 = 7, &
JENDL_40 = 8
! Temperature treatment method
integer, parameter :: &
TEMPERATURE_NEAREST = 1, &
TEMPERATURE_INTERPOLATION = 2, &
TEMPERATURE_MULTIPOLE = 3
! ============================================================================
! TALLY-RELATED CONSTANTS
@ -401,12 +407,6 @@ module constants
integer, parameter :: ERROR_INT = -huge(0)
real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8
! Energy grid methods
integer, parameter :: &
GRID_NUCLIDE = 1, & ! unique energy grid for each nuclide
GRID_MAT_UNION = 2, & ! material union grids with pointers
GRID_LOGARITHM = 3 ! lethargy mapping
! Running modes
integer, parameter :: &
MODE_FIXEDSOURCE = 1, & ! Fixed source mode

514
src/cross_section.F90
View file

@ -1,7 +1,8 @@
module cross_section
use algorithm, only: binary_search
use constants
use energy_grid, only: grid_method, log_spacing
use energy_grid, only: log_spacing
use error, only: fatal_error
use global
use list_header, only: ListElemInt
@ -14,7 +15,6 @@ module cross_section
use particle_header, only: Particle
use random_lcg, only: prn, future_prn, prn_set_stream
use sab_header, only: SAlphaBeta
use search, only: binary_search
implicit none
@ -37,7 +37,6 @@ contains
! 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
! Set all material macroscopic cross sections to zero
material_xs % total = ZERO
@ -49,89 +48,83 @@ contains
! Exit subroutine if material is void
if (p % material == MATERIAL_VOID) return
mat => materials(p % material)
! Find energy index on energy grid
if (grid_method == GRID_MAT_UNION) then
i_grid = find_energy_index(mat, p % E)
else if (grid_method == GRID_LOGARITHM) then
associate (mat => materials(p % material))
! Find energy index on energy grid
i_grid = int(log(p % E/energy_min_neutron)/log_spacing)
end if
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
! Determine if this material has S(a,b) tables
check_sab = (mat % n_sab > 0)
! Initialize position in i_sab_nuclides
j = 1
! Initialize position in i_sab_nuclides
j = 1
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! CHECK FOR S(A,B) TABLE
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! CHECK FOR S(A,B) TABLE
i_sab = 0
i_sab = 0
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
if (check_sab) then
if (i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
if (check_sab) then
if (i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % threshold_inelastic) i_sab = 0
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0
! Increment position in i_sab_nuclides
j = j + 1
! Increment position in i_sab_nuclides
j = j + 1
! Don't check for S(a,b) tables if there are no more left
if (j > mat % n_sab) check_sab = .false.
! Don't check for S(a,b) tables if there are no more left
if (j > mat % n_sab) check_sab = .false.
end if
end if
end if
! ========================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! ========================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! Determine microscopic cross sections for this nuclide
i_nuclide = mat % nuclide(i)
! Determine microscopic cross sections for this nuclide
i_nuclide = mat % nuclide(i)
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, &
i_grid, p % sqrtkT)
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, &
i_grid, p % sqrtkT)
end if
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
end if
! ========================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! ========================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
! Add contributions to material macroscopic total cross section
material_xs % total = material_xs % total + &
atom_density * micro_xs(i_nuclide) % total
! Add contributions to material macroscopic total cross section
material_xs % total = material_xs % total + &
atom_density * micro_xs(i_nuclide) % total
! Add contributions to material macroscopic scattering cross section
material_xs % elastic = material_xs % elastic + &
atom_density * micro_xs(i_nuclide) % elastic
! Add contributions to material macroscopic scattering cross section
material_xs % elastic = material_xs % elastic + &
atom_density * micro_xs(i_nuclide) % elastic
! Add contributions to material macroscopic absorption cross section
material_xs % absorption = material_xs % absorption + &
atom_density * micro_xs(i_nuclide) % absorption
! Add contributions to material macroscopic absorption cross section
material_xs % absorption = material_xs % absorption + &
atom_density * micro_xs(i_nuclide) % absorption
! Add contributions to material macroscopic fission cross section
material_xs % fission = material_xs % fission + &
atom_density * micro_xs(i_nuclide) % fission
! Add contributions to material macroscopic fission cross section
material_xs % fission = material_xs % fission + &
atom_density * micro_xs(i_nuclide) % fission
! 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
end do
! 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
end do
end associate
end subroutine calculate_xs
@ -140,169 +133,161 @@ 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, &
i_log_union, sqrtkT)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT)
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) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
logical :: use_mp ! true if XS can be calculated with windowed multipole
integer :: i_temp ! index for temperature
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
real(8) :: f ! interp factor on nuclide energy grid
real(8) :: kT ! temperature in MeV
real(8) :: sigT, sigA, sigF ! Intermediate multipole variables
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
! Set pointer to nuclide and material
nuc => nuclides(i_nuclide)
mat => materials(i_mat)
! Check to see if there is multipole data present at this energy
use_mp = .false.
if (nuc % mp_present) then
if (E >= nuc % multipole % start_E/1.0e6_8 .and. &
E <= nuc % multipole % end_E/1.0e6_8) then
use_mp = .true.
end if
end if
! Evaluate multipole or interpolate
if (use_mp) then
! Call multipole kernel
call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF)
micro_xs(i_nuclide) % total = sigT
micro_xs(i_nuclide) % absorption = sigA
micro_xs(i_nuclide) % elastic = sigT - sigA
if (nuc % fissionable) then
micro_xs(i_nuclide) % fission = sigF
micro_xs(i_nuclide) % nu_fission = sigF * nuc % nu(E, EMISSION_TOTAL)
associate (nuc => nuclides(i_nuclide))
! Check to see if there is multipole data present at this energy
use_mp = .false.
if (nuc % mp_present) then
if (E >= nuc % multipole % start_E/1.0e6_8 .and. &
E <= nuc % multipole % end_E/1.0e6_8) then
use_mp = .true.
else
! If using multipole data but outside the RRR, pick the nearest
! temperature. Note that there is no tolerance here, so this
! temperature could be very far off!
i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1)
end if
else
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
! If not using multipole data, do a linear search on temperature
kT = sqrtkT**2
do i_temp = 1, size(nuclides(i_nuclide) % kTs)
if (abs(nuclides(i_nuclide) % kTs(i_temp) - kT) < &
K_BOLTZMANN*temperature_tolerance) exit
end do
end if
! Ensure these values are set
! Note, the only time either is used is in one of 4 places:
! 1. physics.F90 - scatter - For inelastic scatter.
! 2. physics.F90 - sample_fission - For partial fissions.
! 3. tally.F90 - score_general - For tallying on MTxxx reactions.
! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes.
! It is worth noting that none of these occur in the resolved
! resonance range, so the value here does not matter.
micro_xs(i_nuclide) % index_grid = 0
micro_xs(i_nuclide) % interp_factor = ZERO
else
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_MAT_UNION)
! Evaluate multipole or interpolate
if (use_mp) then
! Call multipole kernel
call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF)
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
micro_xs(i_nuclide) % total = sigT
micro_xs(i_nuclide) % absorption = sigA
micro_xs(i_nuclide) % elastic = sigT - sigA
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
if (E < nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
if (nuc % fissionable) then
micro_xs(i_nuclide) % fission = sigF
micro_xs(i_nuclide) % nu_fission = sigF * nuc % nu(E, EMISSION_TOTAL)
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
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), &
i_high - i_low + 1, E) + i_low - 1
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
end if
case (GRID_NUCLIDE)
! Perform binary search on the nuclide energy grid in order to determine
! which points to interpolate between
! Ensure these values are set
! Note, the only time either is used is in one of 4 places:
! 1. physics.F90 - scatter - For inelastic scatter.
! 2. physics.F90 - sample_fission - For partial fissions.
! 3. tally.F90 - score_general - For tallying on MTxxx reactions.
! 4. cross_section.F90 - calculate_urr_xs - For unresolved purposes.
! It is worth noting that none of these occur in the resolved
! resonance range, so the value here does not matter.
micro_xs(i_nuclide) % index_temp = i_temp
micro_xs(i_nuclide) % index_grid = 0
micro_xs(i_nuclide) % interp_factor = ZERO
else
associate (grid => nuc % grid(i_temp), xs => nuc % sum_xs(i_temp))
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
if (E <= nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
if (E < grid % energy(1)) then
i_grid = 1
elseif (E > grid % energy(size(grid % energy))) then
i_grid = size(grid % energy) - 1
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
i_low = grid % grid_index(i_log_union)
i_high = grid % grid_index(i_log_union + 1) + 1
! Perform binary search over reduced range
i_grid = binary_search(grid % energy(i_low:i_high), &
i_high - i_low + 1, E) + i_low - 1
end if
! check for rare case where two energy points are the same
if (grid % energy(i_grid) == grid % energy(i_grid + 1)) &
i_grid = i_grid + 1
! calculate interpolation factor
f = (E - grid % energy(i_grid)) / &
(grid % energy(i_grid + 1) - grid % energy(i_grid))
micro_xs(i_nuclide) % index_temp = i_temp
micro_xs(i_nuclide) % index_grid = i_grid
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
+ f * xs % total(i_grid + 1)
! Calculate microscopic nuclide elastic cross section
micro_xs(i_nuclide) % elastic = (ONE - f) * xs % elastic(i_grid) &
+ f * xs % elastic(i_grid + 1)
! Calculate microscopic nuclide absorption cross section
micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( &
i_grid) + f * xs % absorption(i_grid + 1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) &
+ f * xs % fission(i_grid + 1)
! Calculate microscopic nuclide nu-fission cross section
micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( &
i_grid) + f * xs % nu_fission(i_grid + 1)
end if
end associate
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
! if the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then
if (E > nuc % urr_data(i_temp) % energy(1) .and. E < nuc % &
urr_data(i_temp) % energy(nuc % urr_data(i_temp) % n_energy)) then
call calculate_urr_xs(i_nuclide, i_temp, E)
end if
end select
! check for rare case where two energy points are the same
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
! calculate interpolation factor
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
micro_xs(i_nuclide) % index_grid = i_grid
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
+ f * nuc % total(i_grid+1)
! Calculate microscopic nuclide elastic cross section
micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) &
+ f * nuc % elastic(i_grid+1)
! Calculate microscopic nuclide absorption cross section
micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( &
i_grid) + f * nuc % absorption(i_grid+1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) &
+ f * nuc % fission(i_grid+1)
! 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)
end if
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E)
! if the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present) then
if (E > nuc % urr_data % energy(1) .and. &
E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then
call calculate_urr_xs(i_nuclide, E)
end if
end if
micro_xs(i_nuclide) % last_E = E
micro_xs(i_nuclide) % last_index_sab = i_sab
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
micro_xs(i_nuclide) % last_E = E
micro_xs(i_nuclide) % last_index_sab = i_sab
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
end associate
end subroutine calculate_nuclide_xs
@ -312,75 +297,85 @@ contains
! whatever data were taken from the normal Nuclide table.
!===============================================================================
subroutine calculate_sab_xs(i_nuclide, i_sab, E)
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
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
real(8), intent(in) :: sqrtkT ! temperature
integer :: i_grid ! index on S(a,b) energy grid
integer :: i_temp ! temperature index
real(8) :: f ! interp factor on S(a,b) energy grid
real(8) :: inelastic ! S(a,b) inelastic cross section
real(8) :: elastic ! S(a,b) elastic cross section
type(SAlphaBeta), pointer :: sab
real(8) :: kT
! Set flag that S(a,b) treatment should be used for scattering
micro_xs(i_nuclide) % index_sab = i_sab
! Determine temperature for S(a,b) table
kT = sqrtkT**2
do i_temp = 1, size(sab_tables(i_sab) % kTs)
if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < &
K_BOLTZMANN*temperature_tolerance) exit
end do
! Get pointer to S(a,b) table
sab => sab_tables(i_sab)
associate (sab => sab_tables(i_sab) % data(i_temp))
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i_grid = 1
f = ZERO
else
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i_grid)) / &
(sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid))
end if
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i_grid = 1
f = ZERO
else
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i_grid)) / &
(sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid))
end if
! Calculate S(a,b) inelastic scattering cross section
inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + &
f * sab % inelastic_sigma(i_grid + 1)
! Calculate S(a,b) inelastic scattering cross section
inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + &
f * sab % inelastic_sigma(i_grid + 1)
! Check for elastic data
if (E < sab % threshold_elastic) then
! Determine whether elastic scattering is given in the coherent or
! incoherent approximation. For coherent, the cross section is
! represented as P/E whereas for incoherent, it is simply P
! Check for elastic data
if (E < sab % threshold_elastic) then
! Determine whether elastic scattering is given in the coherent or
! incoherent approximation. For coherent, the cross section is
! represented as P/E whereas for incoherent, it is simply P
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
if (E < sab % elastic_e_in(1)) then
! If energy is below that of the lowest Bragg peak, the elastic
! cross section will be zero
elastic = ZERO
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
if (E < sab % elastic_e_in(1)) then
! If energy is below that of the lowest Bragg peak, the elastic
! cross section will be zero
elastic = ZERO
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
elastic = sab % elastic_P(i_grid) / E
end if
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
elastic = sab % elastic_P(i_grid) / E
! Determine index on elastic energy grid
if (E < sab % elastic_e_in(1)) then
i_grid = 1
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
end if
! Get interpolation factor for elastic grid
f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
sab%elastic_e_in(i_grid))
! Calculate S(a,b) elastic scattering cross section
elastic = (ONE - f) * sab % elastic_P(i_grid) + &
f * sab % elastic_P(i_grid + 1)
end if
else
! Determine index on elastic energy grid
if (E < sab % elastic_e_in(1)) then
i_grid = 1
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
end if
! Get interpolation factor for elastic grid
f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
sab%elastic_e_in(i_grid))
! Calculate S(a,b) elastic scattering cross section
elastic = (ONE - f) * sab % elastic_P(i_grid) + &
f * sab % elastic_P(i_grid + 1)
! No elastic data
elastic = ZERO
end if
else
! No elastic data
elastic = ZERO
end if
end associate
! Correct total and elastic cross sections
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - &
@ -390,6 +385,9 @@ contains
! Store S(a,b) elastic cross section for sampling later
micro_xs(i_nuclide) % elastic_sab = elastic
! Save temperature index
micro_xs(i_nuclide) % index_temp_sab = i_temp
end subroutine calculate_sab_xs
!===============================================================================
@ -397,9 +395,9 @@ contains
! from probability tables
!===============================================================================
subroutine calculate_urr_xs(i_nuclide, E)
subroutine calculate_urr_xs(i_nuclide, i_temp, E)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_temp ! temperature index
real(8), intent(in) :: E ! energy
integer :: i_energy ! index for energy
@ -414,7 +412,7 @@ contains
micro_xs(i_nuclide) % use_ptable = .true.
associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data)
associate (nuc => nuclides(i_nuclide), urr => nuclides(i_nuclide) % urr_data(i_temp))
! determine energy table
i_energy = 1
do
@ -433,7 +431,7 @@ contains
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
call prn_set_stream(STREAM_URR_PTABLE)
r = future_prn(int(nuc_zaid_dict % get_key(nuc % zaid), 8))
r = future_prn(int(i_nuclide, 8))
call prn_set_stream(STREAM_TRACKING)
i_low = 1
@ -497,10 +495,10 @@ contains
f = micro_xs(i_nuclide) % interp_factor
! Determine inelastic scattering cross section
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
associate (xs => nuc % reactions(nuc % urr_inelastic) % xs(i_temp))
if (i_energy >= xs % threshold) then
inelastic = (ONE - f) * xs % value(i_energy - xs % threshold + 1) + &
f * xs % value(i_energy - xs % threshold + 2)
end if
end associate
end if

2
src/eigenvalue.F90
View file

@ -4,6 +4,7 @@ module eigenvalue
use message_passing
#endif
use algorithm, only: binary_search
use constants, only: ZERO
use error, only: fatal_error, warning
use global
@ -11,7 +12,6 @@ module eigenvalue
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh
use random_lcg, only: prn, set_particle_seed, advance_prn_seed
use search, only: binary_search
use string, only: to_str
implicit none

2
src/endf_header.F90
View file

@ -2,10 +2,10 @@ module endf_header
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG
use hdf5_interface
use search, only: binary_search
implicit none

2
src/energy_distribution.F90
View file

@ -2,12 +2,12 @@ module energy_distribution
use hdf5
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tabulated1D
use hdf5_interface
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a

176
src/energy_grid.F90
View file

@ -13,64 +13,18 @@ module energy_grid
contains
!===============================================================================
! UNIONIZED_GRID creates a unionized energy grid, for the entire problem or for
! each material, composed of the grids from each nuclide in the entire problem,
! or each material, respectively. Right now, the grid for each nuclide is added
! into a linked list one at a time with an effective insertion sort. Could be
! done with a hash for all energy points and then a quicksort at the end (what
! hash function to use?)
!===============================================================================
subroutine unionized_grid()
integer :: i ! index in nuclides array
integer :: j ! index in materials array
type(ListReal) :: list
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
call write_message("Creating unionized energy grid...", 5)
! add grid points for each nuclide in the material
do j = 1, n_materials
mat => materials(j)
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
call add_grid_points(list, nuc % energy)
end do
! set size of unionized material energy grid
mat % n_grid = list % size()
! create allocated array from linked list
allocate(mat % e_grid(mat % n_grid))
do i = 1, mat % n_grid
mat % e_grid(i) = list % get_item(i)
end do
! delete linked list and dictionary
call list % clear()
end do
! Set pointers to unionized energy grid for each nuclide
call grid_pointers()
end subroutine unionized_grid
!===============================================================================
! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding
! indices on a nuclide energy grid
!===============================================================================
subroutine logarithmic_grid()
integer :: i, j, k ! Loop indices
integer :: t ! temperature index
integer :: M ! Number of equally log-spaced bins
real(8) :: E_max ! Maximum energy in MeV
real(8) :: E_min ! Minimum energy in MeV
real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid
type(Nuclide), pointer :: nuc
! Set minimum/maximum energies
E_max = energy_max_neutron
@ -85,123 +39,29 @@ contains
umesh(:) = [(i*log_spacing, i=0, M)]
do i = 1, n_nuclides_total
! Allocate logarithmic mapping for nuclide
nuc => nuclides(i)
allocate(nuc % grid_index(0:M))
associate (nuc => nuclides(i))
do t = 1, size(nuc % grid)
! Allocate logarithmic mapping for nuclide
allocate(nuc % grid(t) % grid_index(0:M))
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(nuc%energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(nuc % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == nuc % n_grid) then
exit
end if
j = j + 1
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(nuc % grid(t) % energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(nuc % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == size(nuc % grid(t) % energy)) exit
j = j + 1
end do
nuc % grid(t) % grid_index(k) = j
end do
end do
nuc % grid_index(k) = j
end do
end associate
end do
deallocate(umesh)
end subroutine logarithmic_grid
!===============================================================================
! ADD_GRID_POINTS adds energy points from the 'energy' array into a linked list
! of points already stored from previous arrays.
!===============================================================================
subroutine add_grid_points(list, energy)
type(ListReal) :: list
real(8), intent(in) :: energy(:)
integer :: i ! index in energy array
integer :: n ! size of energy array
integer :: current ! current index
real(8) :: E ! actual energy value
i = 1
n = size(energy)
! Set current index to beginning of the list
current = 1
do while (i <= n)
E = energy(i)
! If we've reached the end of the grid energy list, add the remaining
! energy points to the end
if (current > list % size()) then
! Finish remaining energies
do while (i <= n)
call list % append(energy(i))
i = i + 1
end do
exit
end if
if (E < list % get_item(current)) then
! Insert new energy in this position
call list % insert(current, E)
! Advance index in linked list and in new energy grid
i = i + 1
current = current + 1
elseif (E == list % get_item(current)) then
! Found the exact same energy, no need to store duplicates so just
! skip and move to next index
i = i + 1
current = current + 1
else
current = current + 1
end if
end do
end subroutine add_grid_points
!===============================================================================
! GRID_POINTERS creates an array of pointers (ints) for each nuclide to link
! each point on the nuclide energy grid to one on a unionized energy grid
!===============================================================================
subroutine grid_pointers()
integer :: i ! loop index for nuclides
integer :: j ! loop index for nuclide energy grid
integer :: k ! loop index for materials
integer :: index_e ! index on union energy grid
real(8) :: union_energy ! energy on union grid
real(8) :: energy ! energy on nuclide grid
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
do k = 1, n_materials
mat => materials(k)
allocate(mat % nuclide_grid_index(mat % n_nuclides, mat % n_grid))
do i = 1, mat % n_nuclides
nuc => nuclides(mat % nuclide(i))
index_e = 1
energy = nuc % energy(index_e)
do j = 1, mat % n_grid
union_energy = mat % e_grid(j)
if (union_energy >= energy .and. index_e < nuc % n_grid) then
index_e = index_e + 1
energy = nuc % energy(index_e)
end if
mat % nuclide_grid_index(i,j) = index_e - 1
end do
end do
end do
end subroutine grid_pointers
end module energy_grid

52
src/geometry.F90
View file

@ -248,14 +248,15 @@ contains
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
! Set the particle material
! Save previous material and temperature
p % last_material = p % material
if (size(c % material) == 1) then
! Only one material for this cell; assign that one to the particle.
p % material = c % material(1)
else
! Distributed instances of this cell have different materials.
! Determine which instance this is and assign the matching material.
p % last_sqrtkT = p % sqrtkT
! Get distributed offset
if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then
! Distributed instances of this cell have different
! materials/temperatures. Determine which instance this is for
! assigning the matching material/temperature.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
@ -276,37 +277,20 @@ contains
end if
end if
end do
p % material = c % material(offset + 1)
end if
! Set the particle temperature
if (size(c % sqrtkT) == 1) then
! Only one temperature for this cell; assign that one to the particle.
p % sqrtkT = c % sqrtkT(1)
! Save the material
if (size(c % material) > 1) then
p % material = c % material(offset + 1)
else
! Distributed instances of this cell have different temperatures.
! Determine which instance this is and assign the matching temp.
distribcell_index = c % distribcell_index
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type == CELL_FILL) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
end if
end if
end do
p % material = c % material(1)
end if
! Save the temperature
if (size(c % sqrtkT) > 1) then
p % sqrtkT = c % sqrtkT(offset + 1)
else
p % sqrtkT = c % sqrtkT(1)
end if
elseif (c % type == CELL_FILL) then CELL_TYPE

12
src/global.F90
View file

@ -99,13 +99,10 @@ module global
! What to assume for expanding natural elements
integer :: default_expand = ENDF_BVII1
! Whether or not windowed multipole cross sections should be used.
logical :: multipole_active = .false.
! Total amount of nuclide ZAID and dictionary of nuclide ZAID and index --
! this is used when sampling unresolved resonance probability tables
integer(8) :: n_nuc_zaid_total
type(DictIntInt) :: nuc_zaid_dict
! Default temperature and method for choosing temperatures
integer :: temperature_method = TEMPERATURE_NEAREST
real(8) :: temperature_tolerance = 10.0_8
real(8) :: temperature_default = 293.6_8
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
@ -430,7 +427,6 @@ module global
! Various output options
logical :: output_summary = .true.
logical :: output_xs = .false.
logical :: output_tallies = .true.
! ============================================================================

62
src/hdf5_interface.F90
View file

@ -73,6 +73,7 @@ module hdf5_interface
module procedure read_attribute_integer_1D
module procedure read_attribute_integer_2D
module procedure read_attribute_string
module procedure read_attribute_string_1D
end interface read_attribute
interface write_attribute
@ -2425,6 +2426,67 @@ contains
call h5tclose_f(memtype, hdf5_err)
end subroutine read_attribute_string
subroutine read_attribute_string_1D(buffer, obj_id, name)
character(*), target, allocatable, intent(inout) :: buffer(:)
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
integer :: hdf5_err
integer(HID_T) :: space_id
integer(HID_T) :: attr_id
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: maxdims(1)
call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err)
if (allocated(buffer)) then
dims(:) = shape(buffer)
else
call h5aget_space_f(attr_id, space_id, hdf5_err)
call h5sget_simple_extent_dims_f(space_id, dims, maxdims, hdf5_err)
allocate(buffer(dims(1)))
call h5sclose_f(space_id, hdf5_err)
end if
call read_attribute_string_1D_explicit(attr_id, dims, buffer)
call h5aclose_f(attr_id, hdf5_err)
end subroutine read_attribute_string_1D
subroutine read_attribute_string_1D_explicit(attr_id, dims, buffer)
integer(HID_T), intent(in) :: attr_id
integer(HSIZE_T), intent(in) :: dims(1)
character(*), target, intent(inout) :: buffer(dims(1))
integer :: hdf5_err
integer(HID_T) :: filetype
integer(HID_T) :: memtype
integer(SIZE_T) :: size
integer(SIZE_T) :: n
type(c_ptr) :: f_ptr
! Make sure buffer is large enough
call h5aget_type_f(attr_id, filetype, hdf5_err)
call h5tget_size_f(filetype, size, hdf5_err)
if (size > len(buffer(1)) + 1) then
print *, size, len(buffer(1))
call fatal_error("Character buffer is not long enough to &
&read HDF5 string array.")
end if
! Get datatype in memory based on Fortran character
n = len(buffer(1))
call h5tcopy_f(H5T_FORTRAN_S1, memtype, hdf5_err)
call h5tset_size_f(memtype, n, hdf5_err)
! Get pointer to start of string
f_ptr = c_loc(buffer(1)(1:1))
call h5aread_f(attr_id, memtype, f_ptr, hdf5_err)
call h5tclose_f(filetype, hdf5_err)
call h5tclose_f(memtype, hdf5_err)
end subroutine read_attribute_string_1D_explicit
subroutine get_shape(obj_id, dims)
integer(HID_T), intent(in) :: obj_id
integer(HSIZE_T), intent(out) :: dims(:)

76
src/initialize.F90
View file

@ -4,7 +4,7 @@ module initialize
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use energy_grid, only: logarithmic_grid, grid_method
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
@ -17,7 +17,7 @@ module initialize
use material_header, only: Material
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: title, header, print_version, write_message, &
print_usage, write_xs_summary, print_plot
print_usage, print_plot
use random_lcg, only: initialize_prng
use state_point, only: load_state_point
use string, only: to_str, starts_with, ends_with, str_to_int
@ -111,20 +111,8 @@ contains
if (run_mode /= MODE_PLOTTING) then
! Construct information needed for nuclear data
if (run_CE) then
! Set undefined cell temperatures to match the material data.
call lookup_material_temperatures()
! Construct unionized or log energy grid for cross-sections
select case (grid_method)
case (GRID_NUCLIDE)
continue
case (GRID_MAT_UNION)
call time_unionize%start()
call unionized_grid()
call time_unionize%stop()
case (GRID_LOGARITHM)
call logarithmic_grid()
end select
! Construct log energy grid for cross-sections
call logarithmic_grid()
else
! Create material macroscopic data for MGXS
call time_read_xs%start()
@ -158,9 +146,6 @@ contains
else
! Write summary information
if (output_summary) call write_summary()
! Write cross section information
if (output_xs) call write_xs_summary()
end if
end if
@ -1005,57 +990,4 @@ contains
end subroutine allocate_offsets
!===============================================================================
! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to
! find their temperatures from material data.
!===============================================================================
subroutine lookup_material_temperatures()
integer :: i, j, k
real(8) :: min_temp
logical :: warning_given
warning_given = .false.
do i = 1, n_cells
! Ignore non-normal cells and cells with defined temperature.
if (cells(i) % type /= CELL_NORMAL) cycle
if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle
! Set the number of temperatures equal to the number of materials.
deallocate(cells(i) % sqrtkT)
allocate(cells(i) % sqrtkT(size(cells(i) % material)))
! Check each of the cell materials for temperature data.
do j = 1, size(cells(i) % material)
! Arbitrarily set void regions to 0K.
if (cells(i) % material(j) == MATERIAL_VOID) then
cells(i) % sqrtkT(j) = ZERO
cycle
end if
associate (mat => materials(cells(i) % material(j)))
! Find the temperature of the coldest nuclide.
min_temp = nuclides(mat % nuclide(1)) % kT
do k = 2, mat % n_nuclides
! Warn the user if the nuclides don't have identical temperatues.
if (nuclides(mat % nuclide(k)) % kT /= min_temp &
.and. .not. warning_given .and. multipole_active) then
call warning("OpenMC cannot &
&identify the temperature of at least one cell. For the &
&purposes of multipole cross section evaluations, all cells &
&with unknown temperature will be set to the coldest &
&temperature found in the nuclear data for that cell's &
&material")
warning_given = .true.
end if
min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT)
end do
! Set the temperature for this cell instance.
cells(i) % sqrtkT(j) = sqrt(min_temp)
end associate
end do
end do
end subroutine lookup_material_temperatures
end module initialize

291
src/input_xml.F90
View file

@ -2,6 +2,7 @@ module input_xml
use hdf5
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
@ -23,7 +24,8 @@ module input_xml
use set_header, only: SetChar
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with, tokenize, split_string
starts_with, ends_with, tokenize, split_string, &
zero_padded
use tally_header, only: TallyObject
use tally_filter
use tally_initialize, only: add_tallies
@ -359,26 +361,6 @@ contains
! Copy random number seed if specified
if (check_for_node(doc, "seed")) call get_node_value(doc, "seed", seed)
! Energy grid methods
if (check_for_node(doc, "energy_grid")) then
call get_node_value(doc, "energy_grid", temp_str)
else
temp_str = 'logarithm'
end if
select case (trim(temp_str))
case ('nuclide')
grid_method = GRID_NUCLIDE
case ('material-union', 'union')
grid_method = GRID_MAT_UNION
if (trim(temp_str) == 'union') &
call warning('Energy grids will be unionized by material. Global&
& energy grid unionization is no longer an allowed option.')
case ('logarithm', 'logarithmic', 'log')
grid_method = GRID_LOGARITHM
case default
call fatal_error("Unknown energy grid method: " // trim(temp_str))
end select
! Number of bins for logarithmic grid
if (check_for_node(doc, "log_grid_bins")) then
call get_node_value(doc, "log_grid_bins", n_log_bins)
@ -978,14 +960,6 @@ contains
trim(temp_str) == '0') output_summary = .false.
end if
! Check for cross sections option
if (check_for_node(node_output, "cross_sections")) then
call get_node_value(node_output, "cross_sections", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. &
trim(temp_str) == '1') output_xs = .true.
end if
! Check for ASCII tallies output option
if (check_for_node(node_output, "tallies")) then
call get_node_value(node_output, "tallies", temp_str)
@ -1089,20 +1063,6 @@ contains
end select
end if
! Check to see if windowed multipole functionality is requested
if (check_for_node(doc, "use_windowed_multipole")) then
call get_node_value(doc, "use_windowed_multipole", temp_str)
select case (to_lower(temp_str))
case ('true', '1')
multipole_active = .true.
case ('false', '0')
multipole_active = .false.
case default
call fatal_error("Unrecognized value for <use_windowed_multipole> in &
&settings.xml")
end select
end if
call get_node_list(doc, "volume_calc", node_vol_list)
n = get_list_size(node_vol_list)
allocate(volume_calcs(n))
@ -1111,6 +1071,27 @@ contains
call volume_calcs(i) % from_xml(node_vol)
end do
! Get temperature settings
if (check_for_node(doc, "temperature_default")) then
call get_node_value(doc, "temperature_default", temperature_default)
end if
if (check_for_node(doc, "temperature_method")) then
call get_node_value(doc, "temperature_method", temp_str)
select case (to_lower(temp_str))
case ('nearest')
temperature_method = TEMPERATURE_NEAREST
case ('interpolation')
temperature_method = TEMPERATURE_INTERPOLATION
case ('multipole')
temperature_method = TEMPERATURE_MULTIPOLE
case default
call fatal_error("Unknown temperature method: " // trim(temp_str))
end select
end if
if (check_for_node(doc, "temperature_tolerance")) then
call get_node_value(doc, "temperature_tolerance", temperature_tolerance)
end if
! Close settings XML file
call close_xmldoc(doc)
@ -2050,6 +2031,9 @@ contains
integer :: i, j
type(DictCharInt) :: library_dict
type(Library), allocatable :: libraries(:)
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
real(8), allocatable :: material_temps(:)
if (run_CE) then
call read_ce_cross_sections_xml(libraries)
@ -2076,12 +2060,18 @@ contains
end if
! Parse data from materials.xml
call read_materials_xml(libraries, library_dict)
call read_materials_xml(libraries, library_dict, material_temps)
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures(material_temps)
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs%start()
call read_ce_cross_sections(libraries, library_dict)
call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
call time_read_xs%stop()
end if
@ -2092,9 +2082,10 @@ contains
call library_dict % clear()
end subroutine read_materials
subroutine read_materials_xml(libraries, library_dict)
subroutine read_materials_xml(libraries, library_dict, material_temps)
type(Library), intent(in) :: libraries(:)
type(DictCharInt), intent(inout) :: library_dict
real(8), allocatable, intent(out) :: material_temps(:)
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
@ -2110,7 +2101,6 @@ contains
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
character(20) :: name ! name of isotope, e.g. 92235.03c
character(6) :: default_temperature ! Default temperature, e.g., '300K'
character(MAX_WORD_LEN) :: units ! units on density
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
@ -2144,25 +2134,13 @@ contains
! Parse materials.xml file
call open_xmldoc(doc, filename)
! Copy default temperature
if (check_for_node(doc, "default_temperature")) then
call get_node_value(doc, "default_temperature", default_temperature)
else if (.not. run_CE) then
! FIXME This is only necessary while MG mode does not have a
! temperature dependent library implementation.
! Set a default for MG mode to allow MG libraries to not include
! temperatures
default_temperature = '294K'
else
default_temperature = ''
end if
! Get pointer to list of XML <material>
call get_node_list(doc, "material", node_mat_list)
! Allocate cells array
n_materials = get_list_size(node_mat_list)
allocate(materials(n_materials))
allocate(material_temps(n_materials))
! Initialize count for number of nuclides/S(a,b) tables
index_nuclide = 0
@ -2192,14 +2170,11 @@ contains
call get_node_value(node_mat, "name", mat % name)
end if
! Copy material temperature
! Get material default temperature
if (check_for_node(node_mat, "temperature")) then
call get_node_value(node_mat, "temperature", mat % temperature)
else if (default_temperature /= '') then
mat % temperature = default_temperature
call get_node_value(node_mat, "temperature", material_temps(i))
else
call fatal_error("Must specify either a material temperature or a &
&default temperature")
material_temps(i) = ERROR_REAL
end if
! =======================================================================
@ -5719,10 +5694,10 @@ contains
ASSIGN_SAB: do k = 1, size(mat % i_sab_tables)
! In order to know which nuclide the S(a,b) table applies to, we need
! to search through the list of nuclides for one which has a matching
! zaid
! name
associate (sab => sab_tables(mat % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(mat % nuclide)
if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then
if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then
mat % i_sab_nuclides(k) = j
exit FIND_NUCLIDE
end if
@ -5774,16 +5749,16 @@ contains
end do
end subroutine assign_sab_tables
subroutine read_ce_cross_sections(libraries, library_dict)
subroutine read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
type(Library), intent(in) :: libraries(:)
type(DictCharInt), intent(inout) :: library_dict
type(VectorReal), intent(in) :: nuc_temps(:)
type(VectorReal), intent(in) :: sab_temps(:)
integer :: i, j
integer :: i_library
integer :: i_nuclide
integer :: i_sab
integer :: index_nuc_zaid ! index in nuclide ZAID
integer :: zaid ! ZAID of nuclide
integer(HID_T) :: file_id
integer(HID_T) :: group_id
logical :: mp_found ! if windowed multipole libraries were found
@ -5796,8 +5771,6 @@ contains
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
index_nuc_zaid = 0
! Read cross sections
do i = 1, size(materials)
do j = 1, size(materials(i) % names)
@ -5813,8 +5786,8 @@ contains
! Read nuclide data from HDF5
file_id = file_open(libraries(i_library) % path, 'r')
group_id = open_group(file_id, name)
call nuclides(i_nuclide) % from_hdf5(group_id, &
materials(i) % temperature)
call nuclides(i_nuclide) % from_hdf5(group_id, nuc_temps(i_nuclide), &
temperature_method, temperature_tolerance)
call close_group(group_id)
call file_close(file_id)
@ -5824,23 +5797,19 @@ contains
! Determine if minimum/maximum energy for this nuclide is greater/less
! than the previous
energy_min_neutron = max(energy_min_neutron, nuclides(i_nuclide) % energy(1))
energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % energy(&
size(nuclides(i_nuclide) % energy)))
if (size(nuclides(i_nuclide) % grid) >= 1) then
energy_min_neutron = max(energy_min_neutron, &
nuclides(i_nuclide) % grid(1) % energy(1))
energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % &
grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy)))
end if
! Add name and alias to dictionary
call already_read % add(name)
! Construct dictionary mapping nuclide zaids to [1,N] -- used for
! unresolved resonance probability tables
zaid = nuclides(i_nuclide) % zaid
if (.not. nuc_zaid_dict % has_key(zaid)) then
index_nuc_zaid = index_nuc_zaid + 1
call nuc_zaid_dict % add_key(zaid, index_nuc_zaid)
end if
! Read multipole file into the appropriate entry on the nuclides array
if (multipole_active) call read_multipole_data(i_nuclide)
if (temperature_method == TEMPERATURE_MULTIPOLE) &
call read_multipole_data(i_nuclide)
end if
! Check if material is fissionable
@ -5868,8 +5837,8 @@ contains
! Read S(a,b) data from HDF5
file_id = file_open(libraries(i_library) % path, 'r')
group_id = open_group(file_id, name)
call sab_tables(i_sab) % from_hdf5(group_id, &
materials(i) % temperature)
call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), &
temperature_tolerance)
call close_group(group_id)
call file_close(file_id)
@ -5879,14 +5848,13 @@ contains
end do
end do
n_nuc_zaid_total = index_nuc_zaid
! Associate S(a,b) tables with specific nuclides
call assign_sab_tables()
! Show which nuclide results in lowest energy for neutron transport
do i = 1, size(nuclides)
if (nuclides(i) % energy(nuclides(i) % n_grid) == energy_max_neutron) then
if (nuclides(i) % grid(1) % energy(size(nuclides(i) % grid(1) % energy)) &
== energy_max_neutron) then
call write_message("Maximum neutron transport energy: " // &
trim(to_str(energy_max_neutron)) // " MeV for " // &
trim(adjustl(nuclides(i) % name)), 6)
@ -5895,7 +5863,7 @@ contains
end do
! If the user wants multipole, make sure we found a multipole library.
if (multipole_active) then
if (temperature_method == TEMPERATURE_MULTIPOLE) then
mp_found = .false.
do i = 1, size(nuclides)
if (nuclides(i) % mp_present) then
@ -5911,6 +5879,107 @@ contains
end subroutine read_ce_cross_sections
!===============================================================================
! ASSIGN_TEMPERATURES If any cells have undefined temperatures, try to find
! their temperatures from material or global default temperatures
!===============================================================================
subroutine assign_temperatures(material_temps)
real(8), intent(in) :: material_temps(:)
integer :: i, j
integer :: i_material
do i = 1, n_cells
! Ignore non-normal cells and cells with defined temperature.
if (cells(i) % material(1) == NONE) cycle
if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle
! Set the number of temperatures equal to the number of materials.
deallocate(cells(i) % sqrtkT)
allocate(cells(i) % sqrtkT(size(cells(i) % material)))
! Check each of the cell materials for temperature data.
do j = 1, size(cells(i) % material)
! Arbitrarily set void regions to 0K.
if (cells(i) % material(j) == MATERIAL_VOID) then
cells(i) % sqrtkT(j) = ZERO
cycle
end if
! Use material default or global default temperature
i_material = material_dict % get_key(cells(i) % material(j))
if (material_temps(i_material) /= ERROR_REAL) then
cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * &
material_temps(i_material))
else
cells(i) % sqrtkT(j) = sqrt(K_BOLTZMANN * temperature_default)
end if
end do
end do
end subroutine assign_temperatures
!===============================================================================
! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table
! appears at in the model. Later, this list is used to determine the actual
! temperatures to read (which may be different if interpolation is used)
!===============================================================================
subroutine get_temperatures(nuc_temps, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), allocatable, intent(out) :: sab_temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_sab ! index in S(a,b) array
integer :: i_material
real(8) :: temperature ! temperature in Kelvin
allocate(nuc_temps(n_nuclides_total))
allocate(sab_temps(n_sab_tables))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(cells(i) % material(j))
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in nuc_temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
! Add temperature if it hasn't already been added
if (find(nuc_temps(i_nuclide), temperature) == -1) then
call nuc_temps(i_nuclide) % push_back(temperature)
end if
end do NUC_NAMES_LOOP
if (mat % n_sab > 0) then
SAB_NAMES_LOOP: do k = 1, size(mat % sab_names)
! Get index in nuc_temps array
i_sab = sab_dict % get_key(to_lower(mat % sab_names(k)))
! Add temperature if it hasn't already been added
if (find(sab_temps(i_sab), temperature) == -1) then
call sab_temps(i_sab) % push_back(temperature)
end if
end do SAB_NAMES_LOOP
end if
end associate
end do
end do
end subroutine get_temperatures
!===============================================================================
! READ_0K_ELASTIC_SCATTERING
!===============================================================================
@ -5927,6 +5996,9 @@ contains
real(8) :: xs_cdf_sum
character(MAX_WORD_LEN) :: name
type(Nuclide) :: resonant_nuc
type(VectorReal) :: temperature
call temperature % push_back(ZERO)
do i = 1, size(nuclides_0K)
if (nuc % name == nuclides_0K(i) % nuclide) then
@ -5946,13 +6018,14 @@ contains
! Read nuclide data from HDF5
file_id = file_open(libraries(i_library) % path, 'r')
group_id = open_group(file_id, name)
call resonant_nuc % from_hdf5(group_id, '0K')
call resonant_nuc % from_hdf5(group_id, temperature, &
TEMPERATURE_NEAREST, 1000.0_8)
call close_group(group_id)
call file_close(file_id)
! Copy 0K energy grid and elastic scattering cross section
call move_alloc(TO=nuc % energy_0K, FROM=resonant_nuc % energy)
call move_alloc(TO=nuc % elastic_0K, FROM=resonant_nuc % elastic)
call move_alloc(TO=nuc % energy_0K, FROM=resonant_nuc % grid(1) % energy)
call move_alloc(TO=nuc % elastic_0K, FROM=resonant_nuc % sum_xs(1) % elastic)
nuc % n_grid_0K = size(nuc % energy_0K)
! Build CDF for 0K elastic scattering
@ -5987,18 +6060,22 @@ contains
integer, intent(in) :: i_table ! index in nuclides/sab_tables
integer :: i
logical :: file_exists ! Does multipole library exist?
character(7) :: readable ! Is multipole library readable?
character(6) :: zaid_string ! String of the ZAID
character(MAX_FILE_LEN+9) :: filename ! Path to multipole xs library
character(MAX_FILE_LEN) :: filename ! Path to multipole xs library
! For the time being, and I know this is a bit hacky, we just assume
! that the file will be zaid.h5.
! that the file will be ZZZAAAmM.h5.
associate (nuc => nuclides(i_table))
write(zaid_string, '(I6.6)') nuc % zaid
filename = trim(path_multipole) // zaid_string // ".h5"
if (nuc % metastable > 0) then
filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // &
trim(zero_padded(nuc % A, 3)) // 'm' // &
trim(to_str(nuc % metastable)) // ".h5"
else
filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // &
trim(zero_padded(nuc % A, 3)) // ".h5"
end if
! Check if Multipole library exists and is readable
inquire(FILE=filename, EXIST=file_exists, READ=readable)
@ -6019,16 +6096,6 @@ contains
call multipole_read(filename, nuc % multipole, i_table)
nuc % mp_present = .true.
! Recreate nu-fission cross section
if (nuc % fissionable) then
do i = 1, size(nuc % energy)
nuc % nu_fission(i) = nuc % nu(nuc % energy(i), EMISSION_TOTAL) * &
nuc % fission(i)
end do
else
nuc % nu_fission(:) = ZERO
end if
end associate
end subroutine read_multipole_data

3
src/material_header.F90
View file

@ -13,9 +13,6 @@ module material_header
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
character(6) :: temperature ! Temperature of the material
! as presented in the HDF5 library;
! e.g., "300K"
! Energy grid information
integer :: n_grid ! # of union material grid points

10
src/mesh.F90
View file

@ -1,14 +1,14 @@
module mesh
use constants
use global
use mesh_header
use search, only: binary_search
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants
use global
use mesh_header
implicit none
contains

174
src/mgxs_header.F90
View file

@ -19,7 +19,6 @@ module mgxs_header
type, abstract :: Mgxs
character(len=104) :: name ! name of dataset, e.g. 92235.03c
integer :: zaid ! Z and A identifier, e.g. 92235
real(8) :: awr ! Atomic Weight Ratio
real(8) :: kT ! temperature in MeV (k*T)
@ -29,7 +28,6 @@ module mgxs_header
contains
procedure(mgxs_init_file_), deferred :: init_file ! Initialize the data
procedure(mgxs_print_), deferred :: print ! Writes object info
procedure(mgxs_get_xs_), deferred :: get_xs ! Get the requested xs
procedure(mgxs_combine_), deferred :: combine ! initializes object
! Sample the outgoing energy from a fission event
@ -64,12 +62,6 @@ module mgxs_header
integer, intent(in) :: max_order ! Maximum requested order
end subroutine mgxs_init_file_
subroutine mgxs_print_(this, unit)
import Mgxs
class(Mgxs),intent(in) :: this
integer, optional, intent(in) :: unit
end subroutine mgxs_print_
pure function mgxs_get_xs_(this,xstype,gin,gout,uvw,mu) result(xs)
import Mgxs
class(Mgxs), intent(in) :: this
@ -150,7 +142,6 @@ module mgxs_header
contains
procedure :: init_file => mgxsiso_init_file ! Initialize Nuclidic MGXS Data
procedure :: print => mgxsiso_print ! Writes nuclide info
procedure :: get_xs => mgxsiso_get_xs ! Gets Size of Data w/in Object
procedure :: combine => mgxsiso_combine ! inits object
procedure :: sample_fission_energy => mgxsiso_sample_fission_energy
@ -181,7 +172,6 @@ module mgxs_header
contains
procedure :: init_file => mgxsang_init_file ! Initialize Nuclidic MGXS Data
procedure :: print => mgxsang_print ! Writes nuclide info
procedure :: get_xs => mgxsang_get_xs ! Gets Size of Data w/in Object
procedure :: combine => mgxsang_combine ! inits object
procedure :: sample_fission_energy => mgxsang_sample_fission_energy
@ -211,11 +201,6 @@ module mgxs_header
else
this % kT = ZERO
end if
if (check_for_node(node_xsdata, "zaid")) then
call get_node_value(node_xsdata, "zaid", this % zaid)
else
this % zaid = 0
end if
if (check_for_node(node_xsdata, "awr")) then
call get_node_value(node_xsdata, "awr", this % awr)
else
@ -957,164 +942,6 @@ module mgxs_header
end subroutine mgxsang_init_file
!===============================================================================
! MGXS*_PRINT displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine mgxs_print(this, unit_)
class(Mgxs), intent(in) :: this
integer, intent(in) :: unit_
character(MAX_LINE_LEN) :: temp_str
! Basic nuclide information
write(unit_,*) 'MGXS Entry: ' // trim(this % name)
if (this % zaid > 0) then
write(unit_,*) ' ZAID = ' // trim(to_str(this % zaid))
else if (this % zaid < 0) then
write(unit_,*) ' Material id = ' // trim(to_str(-this % zaid))
end if
if (this % awr > ZERO) then
write(unit_,*) ' AWR = ' // trim(to_str(this % awr))
end if
if (this % kT > ZERO) then
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
end if
if (this % scatt_type == ANGLE_LEGENDRE) then
temp_str = "Legendre"
write(unit_,*) ' Scattering Type = ' // trim(temp_str)
select type(this)
type is (MgxsIso)
temp_str = to_str(size(this % scatter % dist(1) % data,dim=1) - 1)
end select
write(unit_,*) ' Scattering Order = ' // trim(temp_str)
else if (this % scatt_type == ANGLE_HISTOGRAM) then
temp_str = "Histogram"
write(unit_,*) ' Scattering Type = ' // trim(temp_str)
select type(this)
type is (MgxsIso)
temp_str = to_str(size(this % scatter % dist(1) % data,dim=1))
end select
write(unit_,*) ' Num. Distribution Bins = ' // trim(temp_str)
else if (this % scatt_type == ANGLE_TABULAR) then
temp_str = "Tabular"
write(unit_,*) ' Scattering Type = ' // trim(temp_str)
select type(this)
type is (MgxsIso)
temp_str = to_str(size(this % scatter % dist(1) % data,dim=1))
end select
write(unit_,*) ' Num. Distribution Points = ' // trim(temp_str)
end if
write(unit_,*) ' Fissionable = ', this % fissionable
end subroutine mgxs_print
subroutine mgxsiso_print(this, unit)
class(MgxsIso), intent(in) :: this
integer, optional, intent(in) :: unit
integer :: unit_ ! unit to write to
integer :: size_total, size_scattmat, size_mgxs
integer :: gin
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write Basic Nuclide Information
call mgxs_print(this, unit_)
! Determine size of mgxs and scattering matrices
size_scattmat = 0
do gin = 1, size(this % scatter % energy)
size_scattmat = size_scattmat + &
2 * size(this % scatter % energy(gin) % data) + &
size(this % scatter % dist(gin) % data)
end do
size_scattmat = size_scattmat + size(this % scatter % scattxs)
size_scattmat = size_scattmat * 8
size_mgxs = size(this % total) + size(this % absorption) + &
size(this % nu_fission) + size(this % k_fission) + &
size(this % fission) + size(this % chi)
size_mgxs = size_mgxs * 8
! Calculate total memory
size_total = size_scattmat + size_mgxs
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes'
write(unit_,*) ' Scattering Matrices = ' // &
trim(to_str(size_scattmat)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine mgxsiso_print
subroutine mgxsang_print(this, unit)
class(MgxsAngle), intent(in) :: this
integer, optional, intent(in) :: unit
integer :: unit_ ! unit to write to
integer :: size_total, size_scattmat, size_mgxs
integer :: ipol, iazi, gin
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Write Basic Nuclide Information
call mgxs_print(this, unit_)
write(unit_,*) ' # of Polar Angles = ' // trim(to_str(this % n_pol))
write(unit_,*) ' # of Azimuthal Angles = ' // trim(to_str(this % n_azi))
! Determine size of mgxs and scattering matrices
size_scattmat = 0
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, size(this % scatter(iazi, ipol) % obj % energy)
size_scattmat = size_scattmat + &
2 * size(this % scatter(iazi, ipol) % obj % energy(gin) % data) + &
size(this % scatter(iazi, ipol) % obj % dist(gin) % data)
end do
size_scattmat = size_scattmat + &
size(this % scatter(iazi, ipol) % obj % scattxs)
end do
end do
size_scattmat = size_scattmat * 8
size_mgxs = size(this % total) + size(this % absorption) + &
size(this % nu_fission) + size(this % k_fission) + &
size(this % fission) + size(this % chi)
size_mgxs = size_mgxs * 8
! Calculate total memory
size_total = size_scattmat + size_mgxs
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_mgxs)) // ' bytes'
write(unit_,*) ' Scattering Matrices = ' // &
trim(to_str(size_scattmat)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine mgxsang_print
!===============================================================================
! MGXS*_GET_XS returns the requested data cross section data
!===============================================================================
@ -1319,7 +1146,6 @@ module mgxs_header
else
this % name = mat % name
end if
this % zaid = -mat % id
this % fissionable = mat % fissionable
this % scatt_type = scatt_type

110
src/multipole.F90
View file

@ -28,13 +28,8 @@ contains
integer(HID_T) :: group_id
! Intermediate loading components
character(len=10) :: version
integer :: NMT
integer :: i, j
integer, allocatable :: MT(:)
logical :: accumulated_fission
character(len=24) :: MT_n ! Takes the form '/nuclide/reactions/MT???'
integer :: is_fissionable
character(len=10) :: version
associate (nuc => nuclides(i_table))
@ -80,111 +75,8 @@ contains
call read_dataset(multipole % curvefit, group_id, "curvefit")
! Delete ACE pointwise data
call read_dataset(nuc % n_grid, group_id, "n_grid")
deallocate(nuc % energy)
deallocate(nuc % total)
deallocate(nuc % elastic)
deallocate(nuc % fission)
deallocate(nuc % nu_fission)
deallocate(nuc % absorption)
allocate(nuc % energy(nuc % n_grid))
allocate(nuc % total(nuc % n_grid))
allocate(nuc % elastic(nuc % n_grid))
allocate(nuc % fission(nuc % n_grid))
allocate(nuc % nu_fission(nuc % n_grid))
allocate(nuc % absorption(nuc % n_grid))
nuc % total(:) = ZERO
nuc % absorption(:) = ZERO
nuc % fission(:) = ZERO
! Read in new energy axis (converting eV to MeV)
call read_dataset(nuc % energy, group_id, "energy_points")
nuc % energy = nuc % energy / 1.0e6_8
! Get count and list of MT tables
call read_dataset(NMT, group_id, "MT_count")
allocate(MT(NMT))
call read_dataset(MT, group_id, "MT_list")
call close_group(group_id)
accumulated_fission = .false.
! Loop over each MT entry and load it into a reaction.
do i = 1, NMT
write(MT_n, '(A, I3.3)') '/nuclide/reactions/MT', MT(i)
group_id = open_group(file_id, MT_n)
! Each MT needs to be treated slightly differently.
select case (MT(i))
case(ELASTIC)
call read_dataset(nuc % elastic, group_id, "MT_sigma")
nuc % total(:) = nuc % total + nuc % elastic
case(N_FISSION)
call read_dataset(nuc % fission, group_id, "MT_sigma")
nuc % total(:) = nuc % total + nuc % fission
nuc % absorption(:) = nuc % absorption + nuc % fission
accumulated_fission = .true.
case default
! Search through all of our secondary reactions
do j = 1, size(nuc % reactions)
if (nuc % reactions(j) % MT == MT(i)) then
! Match found
! Individual Fission components exist, so remove the combined
! fission cross section.
if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
.or. MT(i) == N_3NF) .and. accumulated_fission) then
nuc % total(:) = nuc % total - nuc % fission
nuc % absorption(:) = nuc % absorption - nuc % fission
nuc % fission(:) = ZERO
accumulated_fission = .false.
end if
deallocate(nuc % reactions(j) % sigma)
allocate(nuc % reactions(j) % sigma(nuc % n_grid))
call read_dataset(nuc % reactions(j) % sigma, &
group_id, "MT_sigma")
call read_dataset(nuc % reactions(j) % Q_value, &
group_id, "Q_value")
call read_dataset(nuc % reactions(j) % threshold, &
group_id, "threshold")
nuc % reactions(j) % threshold = 1 ! TODO: reconsider implications.
nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value &
/ 1.0e6_8
! Accumulate total
if (MT(i) /= N_LEVEL .and. MT(i) <= N_DA) then
nuc % total(:) = nuc % total + nuc % reactions(j) % sigma
end if
! Accumulate absorption
if (MT(i) >= N_GAMMA .and. MT(i) <= N_DA) then
nuc % absorption(:) = nuc % absorption &
+ nuc % reactions(j) % sigma
end if
! Accumulate fission (if needed)
if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
.or. MT(i) == N_3NF) ) then
nuc % fission(:) = nuc % fission + nuc % reactions(j) % sigma
nuc % absorption(:) = nuc % absorption &
+ nuc % reactions(j) % sigma
end if
end if
end do
end select
call close_group(group_id)
end do
! Close file
call file_close(file_id)

524
src/nuclide_header.F90
View file

@ -7,6 +7,7 @@ module nuclide_header
h5lget_name_by_idx_f, H5_INDEX_NAME_F, H5_ITER_INC_F
use h5lt, only: h5ltpath_valid_f
use algorithm, only: sort, find
use constants
use dict_header, only: DictIntInt
use endf, only: reaction_name, is_fission, is_disappearance
@ -20,7 +21,7 @@ module nuclide_header
use product_header, only: AngleEnergyContainer
use reaction_header, only: Reaction
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use stl_vector, only: VectorInt
use stl_vector, only: VectorInt, VectorReal
use string
use urr_header, only: UrrData
use xml_interface
@ -32,29 +33,37 @@ module nuclide_header
! for continuous-energy neutron transport.
!===============================================================================
type :: Nuclide
! Nuclide meta-data
character(20) :: name ! name of nuclide, e.g. U235.71c
integer :: zaid ! Z and A identifier, e.g. 92235
integer :: metastable ! metastable state
real(8) :: awr ! Atomic Weight Ratio
real(8) :: kT ! temperature in MeV (k*T)
! Fission information
logical :: fissionable = .false. ! nuclide is fissionable?
! Energy grid information
integer :: n_grid ! # of nuclide grid points
type EnergyGrid
integer, allocatable :: grid_index(:) ! log grid mapping indices
real(8), allocatable :: energy(:) ! energy values corresponding to xs
end type EnergyGrid
! Microscopic cross sections
type SumXS
real(8), allocatable :: total(:) ! total cross section
real(8), allocatable :: elastic(:) ! elastic scattering
real(8), allocatable :: fission(:) ! fission
real(8), allocatable :: nu_fission(:) ! neutron production
real(8), allocatable :: absorption(:) ! absorption (MT > 100)
real(8), allocatable :: heating(:) ! heating
end type SumXS
type :: Nuclide
! Nuclide meta-data
character(20) :: name ! name of nuclide, e.g. U235.71c
integer :: Z ! atomic number
integer :: A ! mass number
integer :: metastable ! metastable state
real(8) :: awr ! Atomic Weight Ratio
real(8), allocatable :: kTs(:) ! temperature in MeV (k*T)
! Fission information
logical :: fissionable = .false. ! nuclide is fissionable?
! Energy grid for each temperature
type(EnergyGrid), allocatable :: grid(:)
! Microscopic cross sections
type(SumXS), allocatable :: sum_xs(:)
! Resonance scattering info
logical :: resonant = .false. ! resonant scatterer?
@ -77,7 +86,7 @@ module nuclide_header
! Unresolved resonance data
logical :: urr_present = .false.
integer :: urr_inelastic
type(UrrData), pointer :: urr_data => null()
type(UrrData), allocatable :: urr_data(:)
! Multipole data
logical :: mp_present = .false.
@ -94,7 +103,6 @@ module nuclide_header
contains
procedure :: clear => nuclide_clear
procedure :: print => nuclide_print
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
@ -131,6 +139,7 @@ module nuclide_header
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide
integer :: index_temp_sab ! temperature index for sab_tables
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Information for URR probability table use
@ -173,23 +182,25 @@ module nuclide_header
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
if (associated(this % urr_data)) deallocate(this % urr_data)
if (associated(this % multipole)) deallocate(this % multipole)
end subroutine nuclide_clear
subroutine nuclide_from_hdf5(this, group_id, temperature)
class(Nuclide), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
character(len=*), intent(in) :: temperature
subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance)
class(Nuclide), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorReal), intent(in) :: temperature ! list of desired temperatures
integer, intent(in) :: method
real(8), intent(in) :: tolerance
integer :: i
integer :: Z
integer :: A
integer :: n
integer :: storage_type
integer :: max_corder
integer :: n_links
integer :: hdf5_err
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: urr_group, nu_group
integer(HID_T) :: energy_group, energy_dset
integer(HID_T) :: kT_group, kT_dset
@ -201,13 +212,14 @@ module nuclide_header
integer(SIZE_T) :: name_len, name_file_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp
character(MAX_FILE_LEN), allocatable :: temperatures(:)
integer, allocatable :: temperatures_integer(:)
character(6) :: my_temperature
integer :: temperature_integer
type(VectorInt) :: MTs
character(MAX_WORD_LEN) :: temp_str
character(MAX_FILE_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
logical :: exists
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
! Get name of nuclide from group
name_len = len(this % name)
@ -216,62 +228,90 @@ module nuclide_header
! Get rid of leading '/'
this % name = trim(this % name(2:))
call read_attribute(Z, group_id, 'Z')
call read_attribute(A, group_id, 'A')
call read_attribute(this % Z, group_id, 'Z')
call read_attribute(this % A, group_id, 'A')
call read_attribute(this % metastable, group_id, 'metastable')
this % zaid = 1000 * Z + A + 400 * this % metastable
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
kT_group = open_group(group_id, 'kTs')
! Before accessing the temperature data, see if the user-provied temperature
! exists. We can find this out by looking at the datasets within kT_group
temperature_integer = &
str_to_int(temperature(1: len_trim(adjustl(temperature)) - 1))
call get_datasets(kT_group, temperatures)
allocate(temperatures_integer(size(temperatures)))
do i = 1, size(temperatures)
temperatures_integer(i) = &
str_to_int(temperatures(i)(1: len_trim(adjustl(temperatures(i))) - 1))
! Determine temperatures available
call get_datasets(kT_group, dset_names)
allocate(temps_available(size(dset_names)))
do i = 1, size(dset_names)
! Read temperature value
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
temps_available(i) = temps_available(i) / K_BOLTZMANN
end do
select case (method)
case (TEMPERATURE_NEAREST)
! Determine actual temperatures to read
TEMP_LOOP: do i = 1, temperature % size()
temp_desired = temperature % data(i)
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
temp_actual = temps_available(i_closest)
if (abs(temp_actual - temp_desired) < tolerance) then
if (find(temps_to_read, nint(temp_actual)) == -1) then
call temps_to_read % push_back(nint(temp_actual))
! Write warning for resonance scattering data if 0K is not available
if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0) then
call warning(trim(this % name) // " does not contain 0K data &
&needed for resonance scattering options selected. Using &
&data at " // trim(to_str(nint(temp_actual))) // " K instead.")
end if
end if
else
call fatal_error("Nuclear data library does not contain cross sections &
&for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
end if
end do TEMP_LOOP
case (TEMPERATURE_INTERPOLATION)
! TODO: Get bounding temperatures
call fatal_error("Temperature interpolation not yet implemented")
case (TEMPERATURE_MULTIPOLE)
! Add first available temperature
call temps_to_read % push_back(nint(temps_available(1)))
end select
! Sort temperatures to read
call sort(temps_to_read)
n_temperature = temps_to_read % size()
allocate(this % kTs(n_temperature))
allocate(this % grid(n_temperature))
do i = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
! Read exact temperature value
call read_dataset(this % kTs(i), kT_group, trim(temp_str))
! Read energy grid
energy_group = open_group(group_id, 'energy')
energy_dset = open_dataset(energy_group, temp_str)
call get_shape(energy_dset, dims)
allocate(this % grid(i) % energy(int(dims(1), 4)))
call read_dataset(this % grid(i) % energy, energy_dset)
call close_dataset(energy_dset)
call close_group(energy_group)
end do
my_temperature = &
temperatures(minloc(abs(temperature_integer - temperatures_integer), &
dim=1))
! Now print a warning if there is no matching temperature and then use the
! closest temperature
if (temperature /= my_temperature) then
if (temperature == '0K') then
call warning(trim(this % name) // " does not contain 0K data &
&needed for the resonance scattering options selected")
else
call warning(trim(this % name) // " does not contain data at a &
&temperature of " // trim(temperature) // "; using the &
&nearest available temperature of " // trim(my_temperature))
end if
end if
kT_dset = open_dataset(kT_group, my_temperature)
call read_dataset(this % kT, kT_dset)
call close_dataset(kT_dset)
call close_group(kT_group)
! Read energy grid
energy_group = open_group(group_id, 'energy')
energy_dset = open_dataset(energy_group, my_temperature)
call get_shape(energy_dset, dims)
this % n_grid = int(dims(1), 4)
allocate(this % energy(this % n_grid))
call read_dataset(this % energy, energy_dset)
call close_dataset(energy_dset)
call close_group(energy_group)
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
call h5gget_info_f(rxs_group, storage_type, n_links, max_corder, hdf5_err)
do j = 0, n_links - 1
call h5lget_name_by_idx_f(rxs_group, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, &
j, temp, hdf5_err, name_len)
if (starts_with(temp, "reaction_")) then
call MTs % push_back(int(str_to_int(temp(10:12))))
j, temp_str, hdf5_err, name_len)
if (starts_with(temp_str, "reaction_")) then
call MTs % push_back(int(str_to_int(temp_str(10:12))))
end if
end do
@ -281,7 +321,7 @@ module nuclide_header
rx_group = open_group(rxs_group, 'reaction_' // trim(&
zero_padded(MTs % data(i), 3)))
call this % reactions(i) % from_hdf5(rx_group, my_temperature)
call this % reactions(i) % from_hdf5(rx_group, temps_to_read)
call close_group(rx_group)
end do
call close_group(rxs_group)
@ -290,32 +330,42 @@ module nuclide_header
call h5ltpath_valid_f(group_id, 'urr', .true., exists, hdf5_err)
if (exists) then
this % urr_present = .true.
allocate(this % urr_data)
urr_group = open_group(group_id, 'urr/' // trim(my_temperature))
call this % urr_data % from_hdf5(urr_group)
allocate(this % urr_data(n_temperature))
do i = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
! Read probability tables for i-th temperature
urr_group = open_group(group_id, 'urr/' // trim(temp_str))
call this % urr_data(i) % from_hdf5(urr_group)
call close_group(urr_group)
! Check for negative values
if (any(this % urr_data(i) % prob < ZERO)) then
call warning("Negative value(s) found on probability table &
&for nuclide " // this % name // " at " // trim(temp_str))
end if
end do
! if the inelastic competition flag indicates that the inelastic cross
! section should be determined from a normal reaction cross section, we
! need to get the index of the reaction
if (this % urr_data % inelastic_flag > 0) then
do i = 1, size(this % reactions)
if (this % reactions(i) % MT == this % urr_data % inelastic_flag) then
this % urr_inelastic = i
if (n_temperature > 0) then
if (this % urr_data(1) % inelastic_flag > 0) then
do i = 1, size(this % reactions)
if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then
this % urr_inelastic = i
end if
end do
! Abort if no corresponding inelastic reaction was found
if (this % urr_inelastic == NONE) then
call fatal_error("Could not find inelastic reaction specified on &
&unresolved resonance probability table.")
end if
end do
! Abort if no corresponding inelastic reaction was found
if (this % urr_inelastic == NONE) then
call fatal_error("Could not find inelastic reaction specified on &
&unresolved resonance probability table.")
end if
end if
! Check for negative values
if (any(this % urr_data % prob < ZERO)) then
call warning("Negative value(s) found on probability table &
&for nuclide " // this % name)
end if
end if
! Check for nu-total
@ -325,8 +375,8 @@ module nuclide_header
! Read total nu data
total_nu = open_dataset(nu_group, 'yield')
call read_attribute(temp, total_nu, 'type')
select case (temp)
call read_attribute(temp_str, total_nu, 'type')
select case (temp_str)
case ('Tabulated1D')
allocate(Tabulated1D :: this % total_nu)
case ('Polynomial')
@ -346,8 +396,8 @@ module nuclide_header
! Check to see if this is polynomial or tabulated data
fer_dset = open_dataset(fer_group, 'q_prompt')
call read_attribute(temp, fer_dset, 'type')
if (temp == 'Polynomial') then
call read_attribute(temp_str, fer_dset, 'type')
if (temp_str == 'Polynomial') then
! Read the prompt Q-value
allocate(Polynomial :: this % fission_q_prompt)
call this % fission_q_prompt % from_hdf5(fer_dset)
@ -358,7 +408,7 @@ module nuclide_header
fer_dset = open_dataset(fer_group, 'q_recoverable')
call this % fission_q_recov % from_hdf5(fer_dset)
call close_dataset(fer_dset)
else if (temp == 'Tabulated1D') then
else if (temp_str == 'Tabulated1D') then
! Read the prompt Q-value
allocate(Tabulated1D :: this % fission_q_prompt)
call this % fission_q_prompt % from_hdf5(fer_dset)
@ -383,108 +433,125 @@ module nuclide_header
subroutine nuclide_create_derived(this)
class(Nuclide), intent(inout) :: this
integer :: i
integer :: j
integer :: k
integer :: i, j, k
integer :: t
integer :: m
integer :: n
integer :: n_grid
integer :: i_fission
type(ListInt) :: MTs
integer :: n_temperature
type(VectorInt) :: MTs
! Allocate and initialize derived cross sections
allocate(this % total(this % n_grid))
allocate(this % elastic(this % n_grid))
allocate(this % fission(this % n_grid))
allocate(this % nu_fission(this % n_grid))
allocate(this % absorption(this % n_grid))
this % total(:) = ZERO
this % elastic(:) = ZERO
this % fission(:) = ZERO
this % nu_fission(:) = ZERO
this % absorption(:) = ZERO
n_temperature = size(this % kTs)
allocate(this % sum_xs(n_temperature))
do i = 1, n_temperature
! Allocate and initialize derived cross sections
n_grid = size(this % grid(i) % energy)
allocate(this % sum_xs(i) % total(n_grid))
allocate(this % sum_xs(i) % elastic(n_grid))
allocate(this % sum_xs(i) % fission(n_grid))
allocate(this % sum_xs(i) % nu_fission(n_grid))
allocate(this % sum_xs(i) % absorption(n_grid))
this % sum_xs(i) % total(:) = ZERO
this % sum_xs(i) % elastic(:) = ZERO
this % sum_xs(i) % fission(:) = ZERO
this % sum_xs(i) % nu_fission(:) = ZERO
this % sum_xs(i) % absorption(:) = ZERO
end do
i_fission = 0
do i = 1, size(this % reactions)
call MTs % append(this % reactions(i) % MT)
call MTs % push_back(this % reactions(i) % MT)
call this % reaction_index % add_key(this % reactions(i) % MT, i)
associate (rx => this % reactions(i))
j = rx % threshold
n = size(rx % sigma)
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle
if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle
! Skip level cross sections if total is available
if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. MTs % contains(N_P)) cycle
if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. MTs % contains(N_D)) cycle
if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. MTs % contains(N_T)) cycle
if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle
if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. MTs % contains(N_A)) cycle
if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. MTs % contains(N_2N)) cycle
if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. find(MTs, N_P) /= -1) cycle
if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. find(MTs, N_D) /= -1) cycle
if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. find(MTs, N_T) /= -1) cycle
if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. find(MTs, N_3HE) /= -1) cycle
if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. find(MTs, N_A) /= -1) cycle
if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle
! Copy elastic
if (rx % MT == ELASTIC) this % elastic(:) = rx % sigma
do t = 1, n_temperature
j = rx % xs(t) % threshold
n = size(rx % xs(t) % value)
! Add contribution to total cross section
this % total(j:j+n-1) = this % total(j:j+n-1) + rx % sigma
! Copy elastic
if (rx % MT == ELASTIC) this % sum_xs(t) % elastic(:) = rx % xs(t) % value
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
this % absorption(j:j+n-1) = this % absorption(j:j+n-1) + rx % sigma
end if
! Add contribution to total cross section
this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + &
rx % xs(t) % value
! Information about fission reactions
if (rx % MT == N_FISSION) then
allocate(this % index_fission(1))
elseif (rx % MT == N_F) then
allocate(this % index_fission(PARTIAL_FISSION_MAX))
this % has_partial_fission = .true.
end if
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
absorption(j:j+n-1) + rx % xs(t) % value
end if
! Add contribution to fission cross section
if (is_fission(rx % MT)) then
this % fissionable = .true.
this % fission(j:j+n-1) = this % fission(j:j+n-1) + rx % sigma
! Also need to add fission cross sections to absorption
this % absorption(j:j+n-1) = this % absorption(j:j+n-1) + rx % sigma
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to this % fission
if (rx % MT == N_FISSION) deallocate(rx % sigma)
! Keep track of this reaction for easy searching later
i_fission = i_fission + 1
this % index_fission(i_fission) = i
this % n_fission = this % n_fission + 1
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, when the angle/energy
! distribution was uncorrelated, no angle was actually sampled. With
! the refactor, an angle is always sampled for an uncorrelated
! distribution even when no angle distribution exists in the ACE file
! (isotropic is assumed). To preserve the RNG stream, we explicitly
! mark fission reactions so that we avoid the angle sampling.
do k = 1, size(rx % products)
if (rx % products(k) % particle == NEUTRON) then
do m = 1, size(rx % products(k) % distribution)
associate (aedist => rx % products(k) % distribution(m) % obj)
select type (aedist)
type is (UncorrelatedAngleEnergy)
aedist % fission = .true.
end select
end associate
end do
! Information about fission reactions
if (t == 1) then
if (rx % MT == N_FISSION) then
allocate(this % index_fission(1))
elseif (rx % MT == N_F) then
allocate(this % index_fission(PARTIAL_FISSION_MAX))
this % has_partial_fission = .true.
end if
end do
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
end if
end associate
end do
end if
! Add contribution to fission cross section
if (is_fission(rx % MT)) then
this % fissionable = .true.
this % sum_xs(t) % fission(j:j+n-1) = this % sum_xs(t) % &
fission(j:j+n-1) + rx % xs(t) % value
! Also need to add fission cross sections to absorption
this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
absorption(j:j+n-1) + rx % xs(t) % value
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to this % fission
if (rx % MT == N_FISSION) deallocate(rx % xs(t) % value)
! Keep track of this reaction for easy searching later
if (t == 1) then
i_fission = i_fission + 1
this % index_fission(i_fission) = i
this % n_fission = this % n_fission + 1
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, when the angle/energy
! distribution was uncorrelated, no angle was actually sampled. With
! the refactor, an angle is always sampled for an uncorrelated
! distribution even when no angle distribution exists in the ACE file
! (isotropic is assumed). To preserve the RNG stream, we explicitly
! mark fission reactions so that we avoid the angle sampling.
do k = 1, size(rx % products)
if (rx % products(k) % particle == NEUTRON) then
do m = 1, size(rx % products(k) % distribution)
associate (aedist => rx % products(k) % distribution(m) % obj)
select type (aedist)
type is (UncorrelatedAngleEnergy)
aedist % fission = .true.
end select
end associate
end do
end if
end do
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
end if
end if ! fission
end do ! temperature
end associate ! rx
end do ! reactions
! Determine number of delayed neutron precursors
if (this % fissionable) then
@ -497,17 +564,16 @@ module nuclide_header
end if
! Calculate nu-fission cross section
if (this % fissionable) then
do i = 1, size(this % energy)
this % nu_fission(i) = this % nu(this % energy(i), EMISSION_TOTAL) * &
this % fission(i)
end do
else
this % nu_fission(:) = ZERO
end if
! Clear MTs set
call MTs % clear()
do t = 1, n_temperature
if (this % fissionable) then
do i = 1, size(this % sum_xs(t) % fission)
this % sum_xs(t) % nu_fission(i) = this % nu(this % grid(t) % energy(i), &
EMISSION_TOTAL) * this % sum_xs(t) % fission(i)
end do
else
this % sum_xs(t) % nu_fission(:) = ZERO
end if
end do
end subroutine nuclide_create_derived
!===============================================================================
@ -575,86 +641,4 @@ module nuclide_header
end function nuclide_nu
!===============================================================================
! NUCLIDE*_PRINT displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine nuclide_print(this, unit)
class(Nuclide), intent(in) :: this
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Initialize totals
size_urr = 0
size_xs = 0
! Basic nuclide information
write(unit_,*) 'Nuclide ' // trim(this % name)
write(unit_,*) ' zaid = ' // trim(to_str(this % zaid))
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid))
write(unit_,*) ' Fissionable = ', this % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(size(this % reactions)))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, size(this % reactions)
associate (rxn => this % reactions(i))
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
rxn % threshold
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * this % n_grid * 4
! Write information about URR probability tables
size_urr = 0
if (this % urr_present) then
associate(urr => this % urr_data)
write(unit_,*) ' Unresolved resonance probability table:'
write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy))
write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob))
write(unit_,*) ' Interpolation = ' // trim(to_str(urr % interp))
write(unit_,*) ' Inelastic flag = ' // trim(to_str(urr % inelastic_flag))
write(unit_,*) ' Absorption flag = ' // trim(to_str(urr % absorption_flag))
write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth
write(unit_,*) ' Min energy = ', trim(to_str(urr % energy(1)))
write(unit_,*) ' Max energy = ', trim(to_str(urr % energy(urr % n_energy)))
! Calculate memory used by probability tables and add to total
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end associate
end if
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
end subroutine nuclide_print
end module nuclide_header

104
src/output.F90
View file

@ -38,43 +38,56 @@ contains
use omp_lib
#endif
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
' d88P" "Y88b 8888b d8888 d88P Y88b', &
' 888 888 88888b.d88888 888 888', &
' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', &
' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', &
' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', &
' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', &
' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', &
'__________________888______________________________________________________', &
' 888', &
' 888'
write(UNIT=OUTPUT_UNIT, FMT='(/23(A/))') &
' %%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' %%%%%%%%%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%%%%%%%%%', &
' ################## %%%%%%%%%%%%%%%%%%%%%%%', &
' ################### %%%%%%%%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%%%%%%', &
' ##################### %%%%%%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%%', &
' ####################### %%%%%%%%%%%%%%%%%', &
' ###################### %%%%%%%%%%%%%%%%%', &
' #################### %%%%%%%%%%%%%%%%%', &
' ################# %%%%%%%%%%%%%%%%%', &
' ############### %%%%%%%%%%%%%%%%', &
' ############ %%%%%%%%%%%%%%%', &
' ######## %%%%%%%%%%%%%%', &
' %%%%%%%%%%%'
! Write version information
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Copyright: 2011-2016 Massachusetts Institute of Technology'
' | The OpenMC Monte Carlo Code'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License: http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') &
' Copyright | 2011-2016 Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License | http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(10X,"Git SHA1 | ",A)') GIT_SHA1
#endif
! Write the date and time
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') &
time_stamp()
write(UNIT=OUTPUT_UNIT, FMT='(9X,"Date/Time | ",A)') time_stamp()
#ifdef MPI
! Write number of processors
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') &
write(UNIT=OUTPUT_UNIT, FMT='(5X,"MPI Processes | ",A)') &
trim(to_str(n_procs))
#endif
#ifdef _OPENMP
! Write number of OpenMP threads
write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') &
write(UNIT=OUTPUT_UNIT, FMT='(4X,"OpenMP Threads | ",A)') &
trim(to_str(omp_get_max_threads()))
#endif
@ -317,57 +330,6 @@ contains
end subroutine print_particle
!===============================================================================
! WRITE_XS_SUMMARY writes information about each nuclide and S(a,b) table to a
! file called cross_sections.out. This file shows the list of reactions as well
! as information about their secondary angle/energy distributions, how much
! memory is consumed, thresholds, etc.
!===============================================================================
subroutine write_xs_summary()
integer :: i ! loop index
integer :: unit_xs ! cross_sections.out file unit
character(MAX_FILE_LEN) :: path ! path of summary file
! Create filename for log file
path = trim(path_output) // "cross_sections.out"
! Open log file for writing
open(NEWUNIT=unit_xs, FILE=path, STATUS='replace', ACTION='write')
if (run_CE) then
! Write header
call header("CROSS SECTION TABLES", unit=unit_xs)
NUCLIDE_LOOP: do i = 1, n_nuclides_total
! Print information about nuclide
call nuclides(i) % print(unit=unit_xs)
end do NUCLIDE_LOOP
SAB_TABLES_LOOP: do i = 1, n_sab_tables
! Print information about S(a,b) table
call sab_tables(i) % print(unit=unit_xs)
end do SAB_TABLES_LOOP
else
! Write header
call header("MGXS LIBRARY TABLES", unit=unit_xs)
NuclideMG_LOOP: do i = 1, n_nuclides_total
! Print information about nuclide
call nuclides_mg(i) % obj % print(unit=unit_xs)
end do NuclideMG_LOOP
call header("MATERIAL MGXS TABLES", unit=unit_xs)
MATERIAL_LOOP: do i = 1, n_materials
! Print information about Materials
call macro_xs(i) % obj % print(unit=unit_xs)
end do MATERIAL_LOOP
end if
! Close cross section summary file
close(unit_xs)
end subroutine write_xs_summary
!===============================================================================
! PRINT_COLUMNS displays a header listing what physical values will displayed
! below them

2
src/particle_header.F90
View file

@ -88,6 +88,7 @@ module particle_header
! Temperature of the current cell
real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV
real(8) :: last_sqrtKT ! last temperature
! Statistical data
integer :: n_collision ! # of collisions
@ -129,6 +130,7 @@ contains
this % cell_born = NONE
this % material = NONE
this % last_material = NONE
this % last_sqrtkT = NONE
this % wgt = ONE
this % last_wgt = ONE
this % absorb_wgt = ZERO

456
src/physics.F90
View file

@ -1,5 +1,6 @@
module physics
use algorithm, only: binary_search
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
@ -15,7 +16,6 @@ module physics
use physics_common
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
use reaction_header, only: Reaction
use search, only: binary_search
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str
@ -59,7 +59,7 @@ contains
! Advance URR seed stream 'N' times after energy changes
if (p % E /= p % last_E) then
call prn_set_stream(STREAM_URR_PTABLE)
call advance_prn_seed(n_nuc_zaid_total)
call advance_prn_seed(size(nuclides, kind=8))
call prn_set_stream(STREAM_TRACKING)
endif
@ -200,6 +200,7 @@ contains
integer :: i
integer :: i_grid
integer :: i_temp
real(8) :: f
real(8) :: prob
real(8) :: cutoff
@ -219,6 +220,7 @@ contains
end if
! Get grid index and interpolatoin factor and sample fission cdf
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
cutoff = prn() * micro_xs(i_nuclide) % fission
@ -229,13 +231,13 @@ contains
FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
i_reaction = nuc % index_fission(i)
associate (rxn => nuc % reactions(i_reaction))
associate (xs => nuc % reactions(i_reaction) % xs(i_temp))
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
if (i_grid < xs % 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)))
prob = prob + ((ONE - f) * xs % value(i_grid - xs % threshold + 1) &
+ f*(xs % value(i_grid - xs % threshold + 2)))
end associate
! Create fission bank sites if fission occurs
@ -294,6 +296,7 @@ contains
integer, intent(in) :: i_nuc_mat
integer :: i
integer :: i_temp
integer :: i_grid
real(8) :: f
real(8) :: prob
@ -301,6 +304,7 @@ contains
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
real(8) :: kT ! temperature in MeV
type(Nuclide), pointer :: nuc
! copy incoming direction
@ -308,6 +312,7 @@ contains
! Get pointer to nuclide and grid index/interpolation factor
nuc => nuclides(i_nuclide)
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
@ -328,8 +333,15 @@ contains
p % E, p % coord(1) % uvw, p % mu)
else
! Determine temperature
if (temperature_method == TEMPERATURE_MULTIPOLE) then
kT = p % sqrtkT**2
else
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), &
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
end if
@ -352,22 +364,24 @@ contains
&// trim(nuc % name))
end if
associate (rxn => nuc % reactions(i))
associate (rx => 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
if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF &
.or. rx % MT == N_2NF .or. rx % 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
if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
associate (xs => rx % xs(i_temp))
! if energy is below threshold for this reaction, skip it
if (i_grid < xs % 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)*xs % value(i_grid - xs % threshold + 1) &
+ f*(xs % value(i_grid - xs % threshold + 2)))
end associate
end associate
end do
@ -401,9 +415,10 @@ contains
! target.
!===============================================================================
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt)
subroutine elastic_scatter(i_nuclide, rxn, kT, E, uvw, mu_lab, wgt)
integer, intent(in) :: i_nuclide
type(Reaction), intent(in) :: rxn
real(8), intent(in) :: kT ! temperature in MeV
real(8), intent(inout) :: E
real(8), intent(inout) :: uvw(3)
real(8), intent(out) :: mu_lab
@ -430,7 +445,7 @@ contains
! Sample velocity of target nucleus
if (.not. micro_xs(i_nuclide) % use_ptable) then
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
& micro_xs(i_nuclide) % elastic)
micro_xs(i_nuclide) % elastic, kT)
else
v_t = ZERO
end if
@ -494,6 +509,7 @@ contains
integer :: i ! incoming energy bin
integer :: j ! outgoing energy bin
integer :: k ! outgoing cosine bin
integer :: i_temp ! temperature index
integer :: n_energy_out ! number of outgoing energy bins
real(8) :: f ! interpolation factor
real(8) :: r ! used for skewed sampling & continuous
@ -502,7 +518,6 @@ contains
real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j)
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
real(8) :: prob ! probability for sampling Bragg edge
type(SAlphaBeta), pointer :: sab
! Following are needed only for SAB_SECONDARY_CONT scattering
integer :: l ! sampled incoming E bin (is i or i + 1)
real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i
@ -514,213 +529,216 @@ contains
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: r1 ! RNG for outgoing energy
i_temp = micro_xs(i_nuclide) % index_temp_sab
! Get pointer to S(a,b) table
sab => sab_tables(i_sab)
associate (sab => sab_tables(i_sab) % data(i_temp))
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % elastic_sab / &
micro_xs(i_nuclide) % elastic) then
! elastic scattering
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % elastic_sab / &
micro_xs(i_nuclide) % elastic) then
! elastic scattering
! Get index and interpolation factor for elastic grid
if (E < sab % elastic_e_in(1)) then
i = 1
f = ZERO
else
i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
f = (E - sab%elastic_e_in(i)) / &
(sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
end if
! Select treatment based on elastic mode
if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then
! With this treatment, we interpolate between two discrete cosines
! corresponding to neighboring incoming energies. This is used for
! data derived in the incoherent approximation
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_elastic_mu)
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % elastic_mu(k,i)
mu_i1jk = sab % elastic_mu(k,i+1)
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then
! This treatment is used for data derived in the coherent
! approximation, i.e. for crystalline structures that have Bragg
! edges.
! Sample a Bragg edge between 1 and i
prob = prn() * sab % elastic_P(i+1)
if (prob < sab % elastic_P(1)) then
k = 1
! Get index and interpolation factor for elastic grid
if (E < sab % elastic_e_in(1)) then
i = 1
f = ZERO
else
k = binary_search(sab % elastic_P(1:i+1), i+1, prob)
i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
f = (E - sab%elastic_e_in(i)) / &
(sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
end if
! Characteristic scattering cosine for this Bragg edge
mu = ONE - TWO*sab % elastic_e_in(k) / E
! Select treatment based on elastic mode
if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then
! With this treatment, we interpolate between two discrete cosines
! corresponding to neighboring incoming energies. This is used for
! data derived in the incoherent approximation
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_elastic_mu)
! Outgoing energy is same as incoming energy -- no need to do anything
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % elastic_mu(k,i)
mu_i1jk = sab % elastic_mu(k,i+1)
else
! Perform inelastic calculations
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i = 1
f = ZERO
else
i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i)) / &
(sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
end if
elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then
! This treatment is used for data derived in the coherent
! approximation, i.e. for crystalline structures that have Bragg
! edges.
! Now that we have an incoming energy bin, we need to determine the
! outgoing energy bin. This will depend on the "secondary energy
! mode". If the mode is 0, then the outgoing energy bin is chosen from a
! set of equally-likely bins. If the mode is 1, then the first
! two and last two bins are skewed to have lower probabilities than the
! other bins (0.1 for the first and last bins and 0.4 for the second and
! second to last bins, relative to a normal bin probability of 1).
! Finally, if the mode is 2, then a continuous distribution (with
! accompanying PDF and CDF is utilized)
if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
j = 1 + int(prn() * sab % n_inelastic_e_out)
elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
! Distribution skewed away from edge points
! Determine number of outgoing energy and angle bins
n_energy_out = sab % n_inelastic_e_out
r = prn() * (n_energy_out - 3)
if (r > ONE) then
! equally likely N-4 middle bins
j = int(r) + 2
elseif (r > 0.6_8) then
! second to last bin has relative probability of 0.4
j = n_energy_out - 1
elseif (r > HALF) then
! last bin has relative probability of 0.1
j = n_energy_out
elseif (r > 0.1_8) then
! second bin has relative probability of 0.4
j = 2
! Sample a Bragg edge between 1 and i
prob = prn() * sab % elastic_P(i+1)
if (prob < sab % elastic_P(1)) then
k = 1
else
! first bin has relative probability of 0.1
j = 1
k = binary_search(sab % elastic_P(1:i+1), i+1, prob)
end if
! Characteristic scattering cosine for this Bragg edge
mu = ONE - TWO*sab % elastic_e_in(k) / E
end if
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
E_ij = sab % inelastic_e_out(j,i)
E_i1j = sab % inelastic_e_out(j,i+1)
! Outgoing energy
E = (1 - f)*E_ij + f*E_i1j
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % inelastic_mu(k,j,i)
mu_i1jk = sab % inelastic_mu(k,j,i+1)
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
else if (sab % secondary_mode == SAB_SECONDARY_CONT) then
! Continuous secondary energy - this is to be similar to
! Law 61 interpolation on outgoing energy
! Sample between ith and (i+1)th bin
r = prn()
if (f > r) then
l = i + 1
else
l = i
end if
! Determine endpoints on grid i
n_energy_out = sab % inelastic_data(i) % n_e_out
E_i_1 = sab % inelastic_data(i) % e_out(1)
E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
! Determine endpoints on grid i + 1
n_energy_out = sab % inelastic_data(i + 1) % n_e_out
E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
E_J = E_i_J + f * (E_i1_J - E_i_J)
! Determine outgoing energy bin
! (First reset n_energy_out to the right value)
n_energy_out = sab % inelastic_data(l) % n_e_out
r1 = prn()
c_j = sab % inelastic_data(l) % e_out_cdf(1)
do j = 1, n_energy_out - 1
c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
if (r1 < c_j1) exit
c_j = c_j1
end do
! check to make sure k is <= n_energy_out - 1
j = min(j, n_energy_out - 1)
! Get the data to interpolate between
E_l_j = sab % inelastic_data(l) % e_out(j)
p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
! Next part assumes linear-linear interpolation in standard
E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
! Find secondary energy (variable E)
frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
if (frac == ZERO) then
E = E_l_j + (r1 - c_j) / p_l_j
else
E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
TWO * frac * (r1 - c_j))) - p_l_j) / frac
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
else
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
end if
! Find angular distribution for closest outgoing energy bin
if (r1 - c_j < c_j1 - r1) then
j = j
else
j = j + 1
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Will use mu from the randomly chosen incoming and closest outgoing
! energy bins
mu = sab % inelastic_data(l) % mu(k, j)
! Outgoing energy is same as incoming energy -- no need to do anything
else
call fatal_error("Invalid secondary energy mode on S(a,b) table " &
&// trim(sab % name))
end if ! (inelastic secondary energy treatment)
end if ! (elastic or inelastic)
! Perform inelastic calculations
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i = 1
f = ZERO
else
i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i)) / &
(sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
end if
! Now that we have an incoming energy bin, we need to determine the
! outgoing energy bin. This will depend on the "secondary energy
! mode". If the mode is 0, then the outgoing energy bin is chosen from a
! set of equally-likely bins. If the mode is 1, then the first
! two and last two bins are skewed to have lower probabilities than the
! other bins (0.1 for the first and last bins and 0.4 for the second and
! second to last bins, relative to a normal bin probability of 1).
! Finally, if the mode is 2, then a continuous distribution (with
! accompanying PDF and CDF is utilized)
if ((sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
(sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED)) then
if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
j = 1 + int(prn() * sab % n_inelastic_e_out)
elseif (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED) then
! Distribution skewed away from edge points
! Determine number of outgoing energy and angle bins
n_energy_out = sab % n_inelastic_e_out
r = prn() * (n_energy_out - 3)
if (r > ONE) then
! equally likely N-4 middle bins
j = int(r) + 2
elseif (r > 0.6_8) then
! second to last bin has relative probability of 0.4
j = n_energy_out - 1
elseif (r > HALF) then
! last bin has relative probability of 0.1
j = n_energy_out
elseif (r > 0.1_8) then
! second bin has relative probability of 0.4
j = 2
else
! first bin has relative probability of 0.1
j = 1
end if
end if
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
E_ij = sab % inelastic_e_out(j,i)
E_i1j = sab % inelastic_e_out(j,i+1)
! Outgoing energy
E = (1 - f)*E_ij + f*E_i1j
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % inelastic_mu(k,j,i)
mu_i1jk = sab % inelastic_mu(k,j,i+1)
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
else if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_CONT) then
! Continuous secondary energy - this is to be similar to
! Law 61 interpolation on outgoing energy
! Sample between ith and (i+1)th bin
r = prn()
if (f > r) then
l = i + 1
else
l = i
end if
! Determine endpoints on grid i
n_energy_out = sab % inelastic_data(i) % n_e_out
E_i_1 = sab % inelastic_data(i) % e_out(1)
E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
! Determine endpoints on grid i + 1
n_energy_out = sab % inelastic_data(i + 1) % n_e_out
E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
E_J = E_i_J + f * (E_i1_J - E_i_J)
! Determine outgoing energy bin
! (First reset n_energy_out to the right value)
n_energy_out = sab % inelastic_data(l) % n_e_out
r1 = prn()
c_j = sab % inelastic_data(l) % e_out_cdf(1)
do j = 1, n_energy_out - 1
c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
if (r1 < c_j1) exit
c_j = c_j1
end do
! check to make sure k is <= n_energy_out - 1
j = min(j, n_energy_out - 1)
! Get the data to interpolate between
E_l_j = sab % inelastic_data(l) % e_out(j)
p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
! Next part assumes linear-linear interpolation in standard
E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
! Find secondary energy (variable E)
frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
if (frac == ZERO) then
E = E_l_j + (r1 - c_j) / p_l_j
else
E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
TWO * frac * (r1 - c_j))) - p_l_j) / frac
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
else
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
end if
! Find angular distribution for closest outgoing energy bin
if (r1 - c_j < c_j1 - r1) then
j = j
else
j = j + 1
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Will use mu from the randomly chosen incoming and closest outgoing
! energy bins
mu = sab % inelastic_data(l) % mu(k, j)
else
call fatal_error("Invalid secondary energy mode on S(a,b) table " &
// trim(sab_tables(i_sab) % name))
end if ! (inelastic secondary energy treatment)
end if ! (elastic or inelastic)
end associate
! Because of floating-point roundoff, it may be possible for mu to be
! outside of the range [-1,1). In these cases, we just set mu to exactly
@ -741,19 +759,19 @@ contains
! implemented here.
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff, kT)
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
real(8), intent(in) :: uvw(3) ! direction cosines
real(8), intent(in) :: v_neut(3) ! neutron velocity
real(8), intent(inout) :: wgt ! particle weight
real(8), intent(in) :: xs_eff ! effective elastic xs at temperature T
real(8), intent(in) :: kT ! equilibrium temperature of target in MeV
real(8) :: awr ! target/neutron mass ratio
real(8) :: kT ! equilibrium temperature of target in MeV
real(8) :: E_rel ! trial relative energy
real(8) :: xs_0K ! 0K xs at E_rel
real(8) :: xs_eff ! effective elastic xs at temperature T
real(8) :: wcf ! weight correction factor
real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1)
real(8) :: E_low ! lowest practical relative energy
@ -782,7 +800,6 @@ contains
character(80) :: sampling_scheme ! method of target velocity sampling
kT = nuc % kT
awr = nuc % awr
! check if nuclide is a resonant scatterer
@ -817,12 +834,12 @@ contains
case ('cxs')
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
case ('wcm')
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! adjust weight as prescribed by the weight correction method (wcm)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
@ -874,7 +891,7 @@ contains
do
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw)
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! perform Doppler broadening rejection correction (dbrc)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
@ -986,13 +1003,13 @@ contains
! can be found in FRA-TM-123.
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
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)
real(8), intent(in) :: kT ! equilibrium temperature of target in MeV
real(8) :: kT ! equilibrium temperature of target in MeV
real(8) :: awr ! target/neutron mass ratio
real(8) :: alpha ! probability of sampling f2 over f1
real(8) :: mu ! cosine of angle between neutron and target vel
@ -1004,7 +1021,6 @@ contains
real(8) :: beta_vt_sq ! (beta * speed of target)^2
real(8) :: vt ! speed of target
kT = nuc % kT
awr = nuc % awr
beta_vn = sqrt(awr * E / kT)

36
src/reaction_header.F90
View file

@ -7,6 +7,7 @@ module reaction_header
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with
implicit none
@ -16,12 +17,16 @@ module reaction_header
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type TemperatureXS
integer :: threshold ! Energy grid index of threshold
real(8), allocatable :: value(:) ! Cross section values
end type TemperatureXS
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
real(8), allocatable :: sigma(:) ! Cross section values
type(TemperatureXS), allocatable :: xs(:)
type(ReactionProduct), allocatable :: products(:)
contains
procedure :: from_hdf5 => reaction_from_hdf5
@ -29,10 +34,10 @@ module reaction_header
contains
subroutine reaction_from_hdf5(this, group_id, temperature)
subroutine reaction_from_hdf5(this, group_id, temperatures)
class(Reaction), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
character(6), intent(in) :: temperature
type(VectorInt), intent(in) :: temperatures
integer :: i
integer :: cm
@ -42,11 +47,12 @@ contains
integer :: n_links
integer :: hdf5_err
integer(HID_T) :: pgroup
integer(HID_T) :: xs, xs_group
integer(HID_T) :: xs, temp_group
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: j
character(MAX_WORD_LEN) :: name
character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K'
call read_attribute(this % Q_value, group_id, 'Q_value')
call read_attribute(this % MT, group_id, 'mt')
@ -54,14 +60,18 @@ contains
this % scatter_in_cm = (cm == 1)
! Read cross section and threshold_idx data
xs_group = open_group(group_id, temperature)
xs = open_dataset(xs_group, 'xs')
call read_attribute(this % threshold, xs, 'threshold_idx')
call get_shape(xs, dims)
allocate(this % sigma(dims(1)))
call read_dataset(this % sigma, xs)
call close_dataset(xs)
call close_group(xs_group)
allocate(this % xs(temperatures % size()))
do i = 1, temperatures % size()
temp_str = trim(to_str(temperatures % data(i))) // "K"
temp_group = open_group(group_id, temp_str)
xs = open_dataset(temp_group, 'xs')
call read_attribute(this % xs(i) % threshold, xs, 'threshold_idx')
call get_shape(xs, dims)
allocate(this % xs(i) % value(dims(1)))
call read_dataset(this % xs(i) % value, xs)
call close_dataset(xs)
call close_group(temp_group)
end do
! Determine number of products
call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err)

6
src/relaxng/materials.rnc
View file

@ -4,7 +4,7 @@ element materials {
(element name { xsd:string { maxLength="52" } } |
attribute name { xsd:string { maxLength="52" } })? &
element temperature { xsd:string { maxLength = "6" } }? &
element temperature { xsd:double }? &
element density {
(element value { xsd:double } | attribute value { xsd:double })? &
@ -41,7 +41,5 @@ element materials {
element sab {
(element name { xsd:string } | attribute name { xsd:string })
}*
}+ &
element default_temperature { xsd:string { maxLength = "6" } }?
}+
}

307
src/relaxng/materials.rng
View file

@ -1,197 +1,186 @@
<?xml version="1.0" encoding="UTF-8"?>
<element name="materials" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<interleave>
<oneOrMore>
<element name="material">
<interleave>
<oneOrMore>
<element name="material">
<interleave>
<choice>
<element name="id">
<data type="int"/>
</element>
<attribute name="id">
<data type="int"/>
</attribute>
</choice>
<optional>
<choice>
<element name="id">
<data type="int"/>
<element name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</element>
<attribute name="id">
<data type="int"/>
<attribute name="name">
<data type="string">
<param name="maxLength">52</param>
</data>
</attribute>
</choice>
<optional>
</optional>
<optional>
<element name="temperature">
<data type="double"/>
</element>
</optional>
<element name="density">
<interleave>
<optional>
<choice>
<element name="value">
<data type="double"/>
</element>
<attribute name="value">
<data type="double"/>
</attribute>
</choice>
</optional>
<choice>
<element name="name">
<element name="units">
<data type="string">
<param name="maxLength">52</param>
<param name="maxLength">10</param>
</data>
</element>
<attribute name="name">
<attribute name="units">
<data type="string">
<param name="maxLength">52</param>
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<element name="temperature">
<data type="string">
<param name="maxLength">6</param>
</data>
</element>
</optional>
<element name="density">
</interleave>
</element>
<zeroOrMore>
<element name="nuclide">
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="value">
<data type="double"/>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="value">
<data type="double"/>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<element name="units">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="units">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
<zeroOrMore>
<element name="nuclide">
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="macroscopic">
</zeroOrMore>
<zeroOrMore>
<element name="macroscopic">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="element">
<interleave>
<choice>
<element name="name">
<data type="string"/>
<data type="string">
<param name="maxLength">2</param>
</data>
</element>
<attribute name="name">
<data type="string"/>
<data type="string">
<param name="maxLength">2</param>
</data>
</attribute>
</choice>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="element">
<interleave>
<optional>
<choice>
<element name="name">
<data type="string">
<param name="maxLength">2</param>
</data>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="name">
<data type="string">
<param name="maxLength">2</param>
</data>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
<optional>
<choice>
<element name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</element>
<attribute name="scattering">
<choice>
<value>data</value>
<value>iso-in-lab</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="sab">
</optional>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
<choice>
<element name="ao">
<data type="double"/>
</element>
<attribute name="ao">
<data type="double"/>
</attribute>
</choice>
<choice>
<element name="wo">
<data type="double"/>
</element>
<attribute name="wo">
<data type="double"/>
</attribute>
</choice>
</choice>
</element>
</zeroOrMore>
</interleave>
</element>
</oneOrMore>
<optional>
<element name="default_temperature">
<data type="string">
<param name="maxLength">6</param>
</data>
</element>
</optional>
</interleave>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="sab">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
</element>
</zeroOrMore>
</interleave>
</element>
</oneOrMore>
</element>

6
src/relaxng/settings.rnc
View file

@ -128,6 +128,12 @@ element settings {
element survival_biasing { xsd:boolean }? &
element temperature_default { xsd:double }? &
element temperature_method { xsd:string }? &
element temperature_tolerance { xsd:double }? &
element threads { xsd:positiveInteger }? &
element trace { list { xsd:positiveInteger+ } }? &

15
src/relaxng/settings.rng
View file

@ -565,6 +565,21 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="temperature_default">
<data type="double"/>
</element>
</optional>
<optional>
<element name="temperature_method">
<data type="string"/>
</element>
</optional>
<optional>
<element name="temperature_tolerance">
<data type="double"/>
</element>
</optional>
<optional>
<element name="threads">
<data type="positiveInteger"/>

428
src/sab_header.F90
View file

@ -2,14 +2,17 @@ module sab_header
use, intrinsic :: ISO_FORTRAN_ENV
use algorithm, only: find, sort
use constants
use dict_header, only: DictIntInt
use distribution_univariate, only: Tabular
use error, only: warning
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use h5lt, only: h5ltpath_valid_f, h5iget_name_f
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets
use secondary_correlated, only: CorrelatedAngleEnergy
use stl_vector, only: VectorInt, VectorReal
use string, only: to_str, str_to_int
implicit none
@ -33,13 +36,7 @@ module sab_header
! of light isotopes such as water, graphite, Be, etc
!===============================================================================
type SAlphaBeta
character(100) :: name ! name of table, e.g. lwtr.10t
real(8) :: awr ! weight of nucleus in neutron masses
real(8) :: kT ! temperature in MeV (k*T)
integer :: n_zaid ! Number of valid zaids
integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012
type SabData
! threshold for S(a,b) treatment (usually ~4 eV)
real(8) :: threshold_inelastic
real(8) :: threshold_elastic = ZERO
@ -48,7 +45,6 @@ module sab_header
integer :: n_inelastic_e_in ! # of incoming E for inelastic
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
real(8), allocatable :: inelastic_e_in(:)
real(8), allocatable :: inelastic_sigma(:)
! The following are used only if secondary_mode is 0 or 1
@ -67,101 +63,35 @@ module sab_header
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
real(8), allocatable :: elastic_mu(:,:)
end type SabData
type SAlphaBeta
character(100) :: name ! name of table, e.g. lwtr.10t
real(8) :: awr ! weight of nucleus in neutron masses
real(8), allocatable :: kTs(:) ! temperatures in MeV (k*T)
character(10), allocatable :: nuclides(:) ! List of valid nuclides
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
! cross sections and distributions at each temperature
type(SabData), allocatable :: data(:)
contains
procedure :: print => salphabeta_print
procedure :: from_hdf5 => salphabeta_from_hdf5
end type SAlphaBeta
contains
!===============================================================================
! PRINT_SAB_TABLE displays information about a S(a,b) table containing data
! describing thermal scattering from bound materials such as hydrogen in water.
!===============================================================================
subroutine salphabeta_print(this, unit)
class(SAlphaBeta), intent(in) :: this
integer, intent(in), optional :: unit
integer :: size_sab ! memory used by S(a,b) table
integer :: unit_ ! unit to write to
integer :: i ! Loop counter for parsing through this % zaid
integer :: char_count ! Counter for the number of characters on a line
! set default unit for writing information
if (present(unit)) then
unit_ = unit
else
unit_ = OUTPUT_UNIT
end if
! Basic S(a,b) table information
write(unit_,*) 'S(a,b) Table ' // trim(this % name)
write(unit_,'(A)',advance="no") ' zaids = '
! Initialize the counter based on the above string
char_count = 11
do i = 1, this % n_zaid
! Deal with a line thats too long
if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma)
! End the line
write(unit_,*) ""
! Add 11 leading blanks
write(unit_,'(A)', advance="no") " "
! reset the counter to 11
char_count = 11
end if
if (i < this % n_zaid) then
! Include a comma
write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i))) // ", "
char_count = char_count + len(trim(to_str(this % zaid(i)))) + 2
else
! Don't include a comma, since we are all done
write(unit_,'(A)',advance="no") trim(to_str(this % zaid(i)))
end if
end do
write(unit_,*) "" ! Move to next line
write(unit_,*) ' awr = ' // trim(to_str(this % awr))
write(unit_,*) ' kT = ' // trim(to_str(this % kT))
! Inelastic data
write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // &
trim(to_str(this % n_inelastic_e_in))
write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // &
trim(to_str(this % n_inelastic_e_out))
write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // &
trim(to_str(this % n_inelastic_mu))
write(unit_,*) ' Threshold for Inelastic = ' // &
trim(to_str(this % threshold_inelastic))
! Elastic data
if (this % n_elastic_e_in > 0) then
write(unit_,*) ' # of Incoming Energies (Elastic) = ' // &
trim(to_str(this % n_elastic_e_in))
write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // &
trim(to_str(this % n_elastic_mu))
write(unit_,*) ' Threshold for Elastic = ' // &
trim(to_str(this % threshold_elastic))
end if
! Determine memory used by S(a,b) table and write out
size_sab = 8 * (this % n_inelastic_e_in * (2 + this % n_inelastic_e_out * &
(1 + this % n_inelastic_mu)) + this % n_elastic_e_in * &
(2 + this % n_elastic_mu))
write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes'
! Blank line at end
write(unit_,*)
end subroutine salphabeta_print
subroutine salphabeta_from_hdf5(this, group_id, temperature)
subroutine salphabeta_from_hdf5(this, group_id, temperature, tolerance)
class(SAlphaBeta), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
character(6), intent(in) :: temperature
type(VectorReal), intent(in) :: temperature ! list of temperatures
real(8), intent(in) :: tolerance
integer :: i, j
integer :: n
integer :: t
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer :: hdf5_err
integer(SIZE_T) :: name_len, name_file_len
integer(HID_T) :: T_group
@ -175,10 +105,13 @@ contains
character(20) :: type
logical :: exists
type(CorrelatedAngleEnergy) :: correlated_dist
character(MAX_FILE_LEN), allocatable :: temperatures(:)
integer, allocatable :: temperatures_integer(:)
character(6) :: my_temperature
integer :: temperature_integer
character(MAX_WORD_LEN) :: temp_str
character(MAX_FILE_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
type(VectorInt) :: temps_to_read
! Get name of table from group
name_len = len(this % name)
@ -188,7 +121,7 @@ contains
this % name = trim(this % name(2:))
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
call read_attribute(this % zaid, group_id, 'zaids')
call read_attribute(this % nuclides, group_id, 'nuclides')
call read_attribute(type, group_id, 'secondary_mode')
select case (type)
case ('equal')
@ -198,168 +131,185 @@ contains
case ('continuous')
this % secondary_mode = SAB_SECONDARY_CONT
end select
this % n_zaid = size(this % zaid)
! Read temperatures
kT_group = open_group(group_id, 'kTs')
! Before accessing the temperature data, see if the user-provied temperature
! exists. We can find this out by looking at the datasets within kT_group
temperature_integer = &
str_to_int(temperature(1: len_trim(adjustl(temperature)) - 1))
call get_datasets(kT_group, temperatures)
allocate(temperatures_integer(size(temperatures)))
do i = 1, size(temperatures)
temperatures_integer(i) = &
str_to_int(temperatures(i)(1: len_trim(adjustl(temperatures(i))) - 1))
! Determine temperatures available
call get_datasets(kT_group, dset_names)
allocate(temps_available(size(dset_names)))
do i = 1, size(dset_names)
! Read temperature value
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
temps_available(i) = temps_available(i) / K_BOLTZMANN
end do
my_temperature = &
temperatures(minloc(abs(temperature_integer - temperatures_integer), &
dim=1))
! Now print a warning if there is no matching temperature and then use the
! closest temperature
if (temperature /= my_temperature) then
if (temperature == '0K') then
call warning(trim(this % name) // " does not contain 0K data &
&needed for the resonance scattering options selected")
! Determine actual temperatures to read
TEMP_LOOP: do i = 1, temperature % size()
temp_desired = temperature % data(i)
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
temp_actual = temps_available(i_closest)
if (abs(temp_actual - temp_desired) < tolerance) then
if (find(temps_to_read, nint(temp_actual)) == -1) then
call temps_to_read % push_back(nint(temp_actual))
end if
else
call warning(trim(this % name) // " does not contain data at a &
&temperature of " // trim(temperature) // "; using the &
&nearest available temperature of " // trim(my_temperature))
call fatal_error("Nuclear data library does not contain cross sections &
&for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
end if
end if
end do TEMP_LOOP
kT_dset = open_dataset(kT_group, my_temperature)
call read_dataset(this % kT, kT_dset)
call close_dataset(kT_dset)
call close_group(kT_group)
! TODO: If using interpolation, add a block to add bounding temperatures for
! each
! Open my_temperature group
T_group = open_group(group_id, my_temperature)
! Sort temperatures to read
call sort(temps_to_read)
! Coherent elastic data
call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err)
if (exists) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
call read_attribute(type, dset_id, 'type')
call get_shape(dset_id, dims2)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
n_temperature = temps_to_read % size()
allocate(this % kTs(n_temperature))
allocate(this % data(n_temperature))
! Set cross section data and type
this % n_elastic_e_in = int(dims2(1), 4)
allocate(this % elastic_e_in(this % n_elastic_e_in))
allocate(this % elastic_P(this % n_elastic_e_in))
this % elastic_e_in(:) = temp(:, 1)
this % elastic_P(:) = temp(:, 2)
select case (type)
case ('tab1')
this % elastic_mode = SAB_ELASTIC_DISCRETE
case ('bragg')
this % elastic_mode = SAB_ELASTIC_EXACT
end select
deallocate(temp)
do t = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(t))) // "K"
! Set elastic threshold
this % threshold_elastic = this % elastic_e_in(this % n_elastic_e_in)
! Read exact temperature value
call read_dataset(this % kTs(t), kT_group, temp_str)
! Read angle distribution
if (this % elastic_mode /= SAB_ELASTIC_EXACT) then
dset_id = open_dataset(elastic_group, 'mu_out')
! Open group for temperature i
T_group = open_group(group_id, temp_str)
! Coherent elastic data
call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err)
if (exists) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
call read_attribute(type, dset_id, 'type')
call get_shape(dset_id, dims2)
this % n_elastic_mu = int(dims2(1), 4)
allocate(this % elastic_mu(dims2(1), dims2(2)))
call read_dataset(this % elastic_mu, dset_id)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
end if
call close_group(elastic_group)
end if
! Set cross section data and type
this % data(t) % n_elastic_e_in = int(dims2(1), 4)
allocate(this % data(t) % elastic_e_in(this % data(t) % n_elastic_e_in))
allocate(this % data(t) % elastic_P(this % data(t) % n_elastic_e_in))
this % data(t) % elastic_e_in(:) = temp(:, 1)
this % data(t) % elastic_P(:) = temp(:, 2)
select case (type)
case ('tab1')
this % data(t) % elastic_mode = SAB_ELASTIC_DISCRETE
case ('bragg')
this % data(t) % elastic_mode = SAB_ELASTIC_EXACT
end select
deallocate(temp)
! Inelastic data
call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err)
if (exists) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')
! Read cross section data
dset_id = open_dataset(inelastic_group, 'xs')
call get_shape(dset_id, dims2)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
! Set cross section data
this % n_inelastic_e_in = int(dims2(1), 4)
allocate(this % inelastic_e_in(this % n_inelastic_e_in))
allocate(this % inelastic_sigma(this % n_inelastic_e_in))
this % inelastic_e_in(:) = temp(:, 1)
this % inelastic_sigma(:) = temp(:, 2)
deallocate(temp)
! Set inelastic threshold
this % threshold_inelastic = this % inelastic_e_in(this % n_inelastic_e_in)
if (this % secondary_mode /= SAB_SECONDARY_CONT) then
! Read energy distribution
dset_id = open_dataset(inelastic_group, 'energy_out')
call get_shape(dset_id, dims2)
this % n_inelastic_e_out = int(dims2(1), 4)
allocate(this % inelastic_e_out(dims2(1), dims2(2)))
call read_dataset(this % inelastic_e_out, dset_id)
call close_dataset(dset_id)
! Set elastic threshold
this % data(t) % threshold_elastic = this % data(t) % elastic_e_in(&
this % data(t) % n_elastic_e_in)
! Read angle distribution
dset_id = open_dataset(inelastic_group, 'mu_out')
call get_shape(dset_id, dims3)
this % n_inelastic_mu = int(dims3(1), 4)
allocate(this % inelastic_mu(dims3(1), dims3(2), dims3(3)))
call read_dataset(this % inelastic_mu, dset_id)
call close_dataset(dset_id)
else
! Read correlated angle-energy distribution
call correlated_dist % from_hdf5(inelastic_group)
if (this % data(t) % elastic_mode /= SAB_ELASTIC_EXACT) then
dset_id = open_dataset(elastic_group, 'mu_out')
call get_shape(dset_id, dims2)
this % data(t) % n_elastic_mu = int(dims2(1), 4)
allocate(this % data(t) % elastic_mu(dims2(1), dims2(2)))
call read_dataset(this % data(t) % elastic_mu, dset_id)
call close_dataset(dset_id)
end if
! Convert to S(a,b) native format
n_energy = size(correlated_dist % energy)
allocate(this % inelastic_data(n_energy))
do i = 1, n_energy
associate (edist => correlated_dist % distribution(i))
! Get number of outgoing energies for incoming energy i
n_energy_out = size(edist % e_out)
this % inelastic_data(i) % n_e_out = n_energy_out
allocate(this % inelastic_data(i) % e_out(n_energy_out))
allocate(this % inelastic_data(i) % e_out_pdf(n_energy_out))
allocate(this % inelastic_data(i) % e_out_cdf(n_energy_out))
! Copy outgoing energy distribution
this % inelastic_data(i) % e_out(:) = edist % e_out
this % inelastic_data(i) % e_out_pdf(:) = edist % p
this % inelastic_data(i) % e_out_cdf(:) = edist % c
do j = 1, n_energy_out
select type (adist => edist % angle(j) % obj)
type is (Tabular)
! On first pass, allocate space for angles
if (j == 1) then
n_mu = size(adist % x)
this % n_inelastic_mu = n_mu
allocate(this % inelastic_data(i) % mu(n_mu, n_energy_out))
end if
! Copy outgoing angles
this % inelastic_data(i) % mu(:, j) = adist % x
end select
end do
end associate
end do
call close_group(elastic_group)
end if
call close_group(inelastic_group)
end if
call close_group(T_group)
! Inelastic data
call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err)
if (exists) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')
! Read cross section data
dset_id = open_dataset(inelastic_group, 'xs')
call get_shape(dset_id, dims2)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
! Set cross section data
this % data(t) % n_inelastic_e_in = int(dims2(1), 4)
allocate(this % data(t) % inelastic_e_in(this % data(t) % n_inelastic_e_in))
allocate(this % data(t) % inelastic_sigma(this % data(t) % n_inelastic_e_in))
this % data(t) % inelastic_e_in(:) = temp(:, 1)
this % data(t) % inelastic_sigma(:) = temp(:, 2)
deallocate(temp)
! Set inelastic threshold
this % data(t) % threshold_inelastic = this % data(t) % inelastic_e_in(&
this % data(t) % n_inelastic_e_in)
if (this % secondary_mode /= SAB_SECONDARY_CONT) then
! Read energy distribution
dset_id = open_dataset(inelastic_group, 'energy_out')
call get_shape(dset_id, dims2)
this % data(t) % n_inelastic_e_out = int(dims2(1), 4)
allocate(this % data(t) % inelastic_e_out(dims2(1), dims2(2)))
call read_dataset(this % data(t) % inelastic_e_out, dset_id)
call close_dataset(dset_id)
! Read angle distribution
dset_id = open_dataset(inelastic_group, 'mu_out')
call get_shape(dset_id, dims3)
this % data(t) % n_inelastic_mu = int(dims3(1), 4)
allocate(this % data(t) % inelastic_mu(dims3(1), dims3(2), dims3(3)))
call read_dataset(this % data(t) % inelastic_mu, dset_id)
call close_dataset(dset_id)
else
! Read correlated angle-energy distribution
call correlated_dist % from_hdf5(inelastic_group)
! Convert to S(a,b) native format
n_energy = size(correlated_dist % energy)
allocate(this % data(t) % inelastic_data(n_energy))
do i = 1, n_energy
associate (edist => correlated_dist % distribution(i))
! Get number of outgoing energies for incoming energy i
n_energy_out = size(edist % e_out)
this % data(t) % inelastic_data(i) % n_e_out = n_energy_out
allocate(this % data(t) % inelastic_data(i) % e_out(n_energy_out))
allocate(this % data(t) % inelastic_data(i) % e_out_pdf(n_energy_out))
allocate(this % data(t) % inelastic_data(i) % e_out_cdf(n_energy_out))
! Copy outgoing energy distribution
this % data(t) % inelastic_data(i) % e_out(:) = edist % e_out
this % data(t) % inelastic_data(i) % e_out_pdf(:) = edist % p
this % data(t) % inelastic_data(i) % e_out_cdf(:) = edist % c
do j = 1, n_energy_out
select type (adist => edist % angle(j) % obj)
type is (Tabular)
! On first pass, allocate space for angles
if (j == 1) then
n_mu = size(adist % x)
this % data(t) % n_inelastic_mu = n_mu
allocate(this % data(t) % inelastic_data(i) % mu(&
n_mu, n_energy_out))
end if
! Copy outgoing angles
this % data(t) % inelastic_data(i) % mu(:, j) = adist % x
end select
end do
end associate
end do
end if
call close_group(inelastic_group)
end if
call close_group(T_group)
end do
call close_group(kT_group)
end subroutine salphabeta_from_hdf5
end module sab_header

2
src/scattdata_header.F90
View file

@ -1,10 +1,10 @@
module scattdata_header
use algorithm, only: binary_search
use constants
use error, only: fatal_error
use math
use random_lcg, only: prn
use search, only: binary_search
implicit none

143
src/search.F90
View file

@ -1,143 +0,0 @@
module search
use constants
implicit none
integer, parameter :: MAX_ITERATION = 64
interface binary_search
module procedure binary_search_real, binary_search_int4, binary_search_int8
end interface binary_search
contains
!===============================================================================
! BINARY_SEARCH performs a binary search of an array to find where a specific
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
pure function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
real(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_real
pure function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
integer, intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int4
pure function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
integer(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
array_index = -1
return
end if
n_iteration = 0
do while (R - L > 1)
! Find values at midpoint
array_index = L + (R - L)/2
if (val >= array(array_index)) then
L = array_index
else
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
array_index = -2
return
end if
end do
array_index = L
end function binary_search_int8
end module search

2
src/secondary_correlated.F90
View file

@ -2,13 +2,13 @@ module secondary_correlated
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, HALF, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: get_shape, read_attribute, open_dataset, &
read_dataset, close_dataset
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This

2
src/secondary_kalbach.F90
View file

@ -2,12 +2,12 @@ module secondary_kalbach
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use hdf5_interface, only: read_attribute, read_dataset, open_dataset, &
close_dataset, get_shape
use random_lcg, only: prn
use search, only: binary_search
!===============================================================================
! KalbachMann represents a correlated angle-energy distribution with the angular

13
src/source.F90
View file

@ -1,5 +1,11 @@
module source
use hdf5, only: HID_T
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use bank_header, only: Bank
use constants
use distribution_univariate, only: Discrete
@ -12,17 +18,10 @@ module source
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream
use search, only: binary_search
use string, only: to_str
use math
use state_point, only: read_source_bank, write_source_bank
#ifdef MPI
use message_passing
#endif
use hdf5, only: HID_T
implicit none
contains

14
src/summary.F90
View file

@ -115,36 +115,31 @@ contains
integer :: i
character(12), allocatable :: nucnames(:)
real(8), allocatable :: awrs(:)
integer, allocatable :: zaids(:)
! Write useful data from nuclide objects
nuclide_group = create_group(file_id, "nuclides")
call write_dataset(nuclide_group, "n_nuclides_total", n_nuclides_total)
! Build array of nuclide names, awrs, and zaids
! Build array of nuclide names and awrs
allocate(nucnames(n_nuclides_total))
allocate(awrs(n_nuclides_total))
allocate(zaids(n_nuclides_total))
do i = 1, n_nuclides_total
if (run_CE) then
nucnames(i) = nuclides(i) % name
awrs(i) = nuclides(i) % awr
zaids(i) = nuclides(i) % zaid
else
nucnames(i) = nuclides_MG(i) % obj % name
awrs(i) = nuclides_MG(i) % obj % awr
zaids(i) = nuclides_MG(i) % obj % zaid
end if
end do
! Write nuclide names, awrs and zaids
! Write nuclide names and awrs
call write_dataset(nuclide_group, "names", nucnames)
call write_dataset(nuclide_group, "awrs", awrs)
call write_dataset(nuclide_group, "zaids", zaids)
call close_group(nuclide_group)
deallocate(nucnames, awrs, zaids)
deallocate(nucnames, awrs)
end subroutine write_nuclides
@ -529,9 +524,6 @@ contains
! Write name for this material
call write_dataset(material_group, "name", m % name)
! Write temperature for this material
call write_dataset(material_group, "temperature", m % temperature)
! Write atom density with units
call write_dataset(material_group, "atom_density", m % density)
call write_attribute_string(material_group, "atom_density", "units", &

52
src/tally.F90
View file

@ -1,5 +1,10 @@
module tally
#ifdef MPI
use message_passing
#endif
use algorithm, only: binary_search
use constants
use error, only: fatal_error
use geometry_header
@ -11,15 +16,10 @@ module tally
use mesh_header, only: RegularMesh
use output, only: header
use particle_header, only: LocalCoord, Particle
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyResult
use tally_filter
#ifdef MPI
use message_passing
#endif
implicit none
integer :: position(N_FILTER_TYPES - 3) = 0 ! Tally map positioning array
@ -88,6 +88,7 @@ contains
integer :: l ! loop index for nuclides in material
integer :: m ! loop index for reactions
integer :: q ! loop index for scoring bins
integer :: i_temp ! temperature index
integer :: i_nuc ! index in nuclides array (from material)
integer :: i_energy ! index in nuclide energy grid
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
@ -887,16 +888,18 @@ contains
if (i_nuclide > 0) 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
f = micro_xs(i_nuclide) % interp_factor
if (i_energy >= rxn % threshold) then
score = ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
! Retrieve temperature and energy grid index and interpolation
! factor
i_temp = micro_xs(i_nuclide) % index_temp
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp))
if (i_energy >= xs % threshold) then
score = ((ONE - f) * xs % value(i_energy - &
xs % threshold + 1) + f * xs % value(i_energy - &
xs % threshold + 2)) * atom_density * flux
end if
end associate
end if
@ -911,15 +914,18 @@ contains
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
f = micro_xs(i_nuc) % interp_factor
if (i_energy >= rxn % threshold) then
score = score + ((ONE - f) * rxn % sigma(i_energy - &
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density_ * flux
! Retrieve temperature and energy grid index and interpolation
! factor
i_temp = micro_xs(i_nuc) % index_temp
i_energy = micro_xs(i_nuc) % index_grid
f = micro_xs(i_nuc) % interp_factor
associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp))
if (i_energy >= xs % threshold) then
score = score + ((ONE - f) * xs % value(i_energy - &
xs % threshold + 1) + f * xs % value(i_energy - &
xs % threshold + 2)) * atom_density_ * flux
end if
end associate
end if

2
src/tally_filter.F90
View file

@ -1,5 +1,6 @@
module tally_filter
use algorithm, only: binary_search
use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION
use dict_header, only: DictIntInt
use geometry_header, only: BASE_UNIVERSE, RectLattice, HexLattice
@ -10,7 +11,6 @@ module tally_filter
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_2d, mesh_intersects_3d
use particle_header, only: Particle
use search, only: binary_search
use string, only: to_str
use tally_filter_header, only: TallyFilter, TallyFilterContainer

7
src/tracking.F90
View file

@ -84,9 +84,10 @@ contains
! Calculate microscopic and macroscopic cross sections
if (run_CE) then
! If the material is the same as the last material and the energy of the
! particle hasn't changed, we don't need to lookup cross sections again.
if (p % material /= p % last_material) call calculate_xs(p)
! If the material is the same as the last material and the temperature
! hasn't changed, we don't need to lookup cross sections again.
if (p % material /= p % last_material .or. &
p % sqrtkT /= p % last_sqrtkT) call calculate_xs(p)
else
! Since the MGXS can be angle dependent, this needs to be done
! After every collision for the MGXS mode

4
tests/test_cmfd_feed/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U235" wo="0.21" />

4
tests/test_cmfd_nofeed/materials.xml
View file

@ -1,9 +1,7 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<!-- Definition of materials -->
<material id="1">
<density value="19" units="g/cc" />
<nuclide name="U235" wo="0.21" />

2
tests/test_complex_cell/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_confidence_intervals/materials.xml
View file

@ -2,7 +2,7 @@
<materials>
<material id="1">
<temperature>294K</temperature>
<temperature>294</temperature>
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>

2
tests/test_density/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="0.1" units="atom/b-cm" />
<nuclide name="U235" ao="1.0" />

2
tests/test_eigenvalue_genperbatch/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_eigenvalue_no_inactive/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_energy_grid/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_energy_grid/settings.xml
View file

@ -1,7 +1,7 @@
<?xml version="1.0"?>
<settings>
<energy_grid>nuclide</energy_grid>
<log_grid_bins>20000</log_grid_bins>
<eigenvalue>
<batches>10</batches>

1
tests/test_energy_laws/materials.xml
View file

@ -1,6 +1,5 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="20" units="g/cc" />
<nuclide name="U233" ao="1.0" />

2
tests/test_entropy/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_filter_distribcell/case-1/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />

3
tests/test_filter_distribcell/case-2/materials.xml
View file

@ -1,9 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Definition of materials -->
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />

2
tests/test_filter_distribcell/case-3/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_filter_distribcell/case-4/materials.xml
View file

@ -1,7 +1,5 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1" name="fuel">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />

2
tests/test_filter_mesh_2d/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_filter_mesh_3d/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<!-- Fuel composition -->
<material id="1">
<density value="10.062" units="g/cm3" />

2
tests/test_fixed_source/materials.xml
View file

@ -1,8 +1,6 @@
<?xml version="1.0"?>
<materials>
<default_temperature>294K</default_temperature>
<material id="1">
<density value="7.5" units="g/cc" />
<nuclide name="O16" ao="1.0" />

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