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

Merge pull request #494 from nelsonag/mg

Browse source 
Addition of MultiGroup Capability
This commit is contained in:
Paul Romano 2016-02-24 06:50:51 -06:00
commit 446262f709
88 changed files with 20739 additions and 1040 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.travis.yml
View file

@ -44,7 +44,7 @@ before_script:
- git clone --branch=master git://github.com/bhermanmit/nndc_xs nndc_xs
- cat nndc_xs/nndc.tar.gza* | tar xzvf -
- rm -rf nndc_xs
- export CROSS_SECTIONS=$PWD/nndc/cross_sections.xml
- export OPENMC_CROSS_SECTIONS=$PWD/nndc/cross_sections.xml
- cd ..
script:

7
docs/source/index.rst
View file

@ -4,9 +4,10 @@ The OpenMC Monte Carlo Code
OpenMC is a Monte Carlo particle transport simulation code focused on neutron
criticality calculations. It is capable of simulating 3D models based on
constructive solid geometry with second-order surfaces. The particle interaction
data is based on ACE format cross sections, also used in the MCNP and Serpent
Monte Carlo codes.
constructive solid geometry with second-order surfaces. OpenMC supports either
continuous-energy or multi-group transport. The continuous-energy
particle interaction data is based on ACE format cross sections, also used
in the MCNP and Serpent Monte Carlo codes.
OpenMC was originally developed by members of the `Computational Reactor Physics
Group`_ at the `Massachusetts Institute of Technology`_ starting

2
docs/source/usersguide/beginners.rst
View file

@ -12,7 +12,7 @@ In a nutshell, OpenMC simulates neutrons moving around randomly in a `nuclear
reactor`_ (or other fissile system). This is what's known as `Monte Carlo`_
simulation. Neutrons are important in nuclear reactors because they are the
particles that induce `fission`_ in uranium and other nuclides. Knowing the
behavior of neutrons allows you to figure out how often and where fission
behavior of neutrons allows you to determine how often and where fission
occurs. The amount of energy released is then directly proportional to the
fission reaction rate since most heat is produced by fission. By simulating many
neutrons (millions or billions), it is possible to determine the average

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

@ -14,6 +14,7 @@ essential aspects of using OpenMC to perform simulations.
beginners
install
input
mgxs_library
output/index
processing
troubleshoot

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

@ -112,9 +112,11 @@ standard deviation.
The ``<cross_sections>`` element has no attributes and simply indicates the path
to an XML cross section listing file (usually named cross_sections.xml). If this
element is absent from the settings.xml file, the :envvar:`CROSS_SECTIONS`
environment variable will be used to find the path to the XML cross section
listing.
element is absent from the settings.xml file, the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable will be used to find the
path to the XML cross section listing when in continuous-energy mode, and the
:envvar:`OPENMC_MG_CROSS_SECTIONS` environment variable will be used in
multi-group mode.
``<cutoff>`` Element
--------------------
@ -212,8 +214,21 @@ cross section values between.
*Default*: logarithm
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
.. _LA-UR-14-24530: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
.. _energy_mode:
``<energy_mode>`` Element
-------------------------
The ``<energy_mode>`` element tells OpenMC if the run-mode should be
continuous-energy or multi-group. Options for entry are: ``continuous-energy``
or ``multi-group``.
*Default*: continuous-energy
``<entropy>`` Element
---------------------
@ -264,6 +279,20 @@ based on the recommended value in LA-UR-14-24530_.
*Default*: 8000
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
``<max_order>`` Element
---------------------------
The ``<max_order>`` element allows the user to set a maximum scattering order
to apply to every nuclide/material in the problem. That is, if the data
library has :math:`P_3` data available, but ``<max_order>`` was set to ``1``,
then, OpenMC will only use up to the :math:`P_1` data.
*Default*: Use the maximum order in the data library
.. note:: This element is not used in the continuous-energy :ref:`energy_mode`.
.. _natural_elements:
``<natural_elements>`` Element
@ -343,6 +372,8 @@ or sub-elements and can be set to either "false" or "true".
*Default*: true
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
``<resonance_scattering>`` Element
----------------------------------
@ -402,6 +433,8 @@ attributes or sub-elements:
*Defaults*: None (scatterer), ARES (method), 0.01 eV (E_min), 1.0 keV (E_max)
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
``<run_cmfd>`` Element
----------------------
@ -579,6 +612,8 @@ variable and whose sub-elements/attributes are as follows:
number :math:`a` that parameterizes the distribution :math:`p(x) dx = c x
e^{-x/a} dx`.
.. note:: The above format should be used even when using the multi-group
:ref:`energy_mode`.
:interpolation:
For a "tabular" distribution, ``interpolation`` can be set to "histogram" or
"linear-linear" thereby specifying how tabular points are to be interpolated.
@ -1184,10 +1219,18 @@ Each ``material`` element can have the following attributes or sub-elements:
``<element>`` sub-elements are to be interpreted as nuclide/element
densities in atom/b-cm, and the total density of the material is taken as
the sum of all nuclides/elements. The "sum" option cannot be used in
conjunction with weight percents.
conjunction with weight percents. The "macro" unit is used with
a ``macroscopic`` quantity to indicate that the density is already included
in the library and thus not needed here. However, if a value is provided
for the ``value``, then this is treated as a number density multiplier on
the macroscopic cross sections in the multi-group data. This can be used,
for example, when perturbing the density slightly.
*Default*: None
.. note:: A ``macroscopic`` quantity can not be used in conjunction with a
``nuclide``, ``element``, or ``sab`` quantity.
:nuclide:
An element with attributes/sub-elements called ``name``, ``xs``, and ``ao``
or ``wo``. The ``name`` attribute is the name of the cross-section for a
@ -1258,6 +1301,24 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: None
:macroscopic:
The ``macroscopic`` element is similar to the ``nuclide`` element, but,
recognizes that some multi-group libraries may be providing material
specific macroscopic cross sections instead of always providing nuclide
specific data like in the continuous-energy case. To that end, the
macroscopic element has attributes/sub-elements called ``name``, and ``xs``.
The ``name`` attribute is the name of the cross-section for a
desired nuclide while the ``xs`` attribute is the cross-section
identifier. One example would be as follows:
.. code-block:: xml
<macroscopic name="UO2" xs="71c" />
.. note:: This element is only used in the multi-group :ref:`energy_mode`.
*Default*: None
.. _IUPAC Isotopic Compositions of the Elements 2009:
http://pac.iupac.org/publications/pac/pdf/2011/pdf/8302x0397.pdf
@ -1344,7 +1405,8 @@ The ``<tally>`` element accepts the following sub-elements:
A list of universes for which the tally should be accumulated.
:energy:
A monotonically increasing list of bounding **pre-collision** energies
In continuous-energy mode, this filter should be provided as a
monotonically increasing list of bounding **pre-collision** energies
for a number of groups. For example, if this filter is specified as
.. code-block:: xml
@ -1354,17 +1416,24 @@ The ``<tally>`` element accepts the following sub-elements:
then two energy bins will be created, one with energies between 0 and
1 MeV and the other with energies between 1 and 20 MeV.
In multi-group mode the bins provided must match group edges
defined in the multi-group library.
:energyout:
A monotonically increasing list of bounding **post-collision**
energies for a number of groups. For example, if this filter is
specified as
In continuous-energy mode, this filter should be provided as a
monotonically increasing list of bounding **post-collision** energies
for a number of groups. For example, if this filter is specified as
.. code-block:: xml
<filter type="energyout" bins="0.0 1.0 20.0" />
then two post-collision energy bins will be created, one with energies
between 0 and 1 MeV and the other with energies between 1 and 20 MeV.
then two post-collision energy bins will be created, one with
energies between 0 and 1 MeV and the other with energies between
1 and 20 MeV.
In multi-group mode the bins provided must match group edges
defined in the multi-group library.
:mu:
A monotonically increasing list of bounding **post-collision** cosines
@ -1448,6 +1517,8 @@ The ``<tally>`` element accepts the following sub-elements:
<filter type="delayedgroup" bins="1 2 3 4 5 6" />
.. note:: This filter type is not used in the multi-group :ref:`energy_mode`.
:nuclides:
If specified, the scores listed will be for particular nuclides, not the
summation of reactions from all nuclides. The format for nuclides should be
@ -1627,7 +1698,8 @@ The ``<tally>`` element accepts the following sub-elements:
|Score | Description |
+======================+===================================================+
|delayed-nu-fission |Total production of delayed neutrons due to |
| |fission. |
| |fission. This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
@ -1656,6 +1728,8 @@ The ``<tally>`` element accepts the following sub-elements:
+----------------------+---------------------------------------------------+
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
| |This score type is not used in the |
| |multi-group :ref:`energy_mode`. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |

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

@ -356,7 +356,7 @@ Testing Build
-------------
If you have ENDF/B-VII.1 cross sections from NNDC_ you can test your build.
Make sure the **CROSS_SECTIONS** environmental variable is set to the
Make sure the **OPENMC_CROSS_SECTIONS** environmental variable is set to the
*cross_sections.xml* file in the *data/nndc* directory.
There are two ways to run tests. The first is to use the Makefile present in
the source directory and run the following:
@ -381,11 +381,17 @@ Cross Section Configuration
---------------------------
In order to run a simulation with OpenMC, you will need cross section data for
each nuclide in your problem. Since OpenMC uses ACE format cross sections, you
can use nuclear data that was processed with NJOY_, such as that distributed
with MCNP_ or Serpent_. Several sources provide free processed ACE data as
described below. The TALYS-based evaluated nuclear data library, TENDL_, is also
openly available in ACE format.
each nuclide or material in your problem. OpenMC can be run in
continuous-energy or multi-group mode.
In continuous-energy mode OpenMC uses ACE format cross sections; in this case
you can use nuclear data that was processed with NJOY_, such as that
distributed with MCNP_ or Serpent_. Several sources provide free processed
ACE data as described below. The TALYS-based evaluated nuclear data library,
TENDL_, is also openly available in ACE format.
In multi-group mode, OpenMC utilizes an XML-based library format which can be
used to describe nuclide- or material-specific quantities.
Using ENDF/B-VII.1 Cross Sections from NNDC
-------------------------------------------
@ -400,9 +406,10 @@ extract, and set up a confiuration file:
cd openmc/data
python get_nndc_data.py
At this point, you should set the :envvar:`CROSS_SECTIONS` environment variable
to the absolute path of the file ``openmc/data/nndc/cross_sections.xml``. This
cross section set is used by the test suite.
At this point, you should set the :envvar:`OPENMC_CROSS_SECTIONS` environment
variable to the absolute path of the file
``openmc/data/nndc/cross_sections.xml``. This cross section set is used by the
test suite.
Using JEFF Cross Sections from OECD/NEA
---------------------------------------
@ -428,8 +435,8 @@ the following steps must be taken:
4. Additionally, you may need to change any occurrences of upper-case "ACE"
within the ``cross_sections.xml`` file to lower-case.
5. Either set the :ref:`cross_sections` in a settings.xml file or the
:envvar:`CROSS_SECTIONS` environment variable to the absolute path of the
``cross_sections.xml`` file.
:envvar:`OPENMC_CROSS_SECTIONS` environment variable to the absolute path of
the ``cross_sections.xml`` file.
Using Cross Sections from MCNP
------------------------------
@ -437,8 +444,9 @@ Using Cross Sections from MCNP
To use cross sections distributed with MCNP, change the <directory> element in
the ``cross_sections.xml`` file in the root directory of the OpenMC distribution
to the location of the MCNP cross sections. Then, either set the
:ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS`
environment variable to the absolute path of the ``cross_sections.xml`` file.
:ref:`cross_sections` in a settings.xml file or the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable to the absolute path of
the ``cross_sections.xml`` file.
Using Cross Sections from Serpent
---------------------------------
@ -446,10 +454,21 @@ Using Cross Sections from Serpent
To use cross sections distributed with Serpent, change the <directory> element
in the ``cross_sections_serpent.xml`` file in the root directory of the OpenMC
distribution to the location of the Serpent cross sections. Then, either set the
:ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS`
environment variable to the absolute path of the ``cross_sections_serpent.xml``
:ref:`cross_sections` in a settings.xml file or the
:envvar:`OPENMC_CROSS_SECTIONS` environment variable to the absolute path of
the ``cross_sections_serpent.xml``
file.
Using Multi-Group Cross Sections
--------------------------------
Multi-group cross section libraries are generally tailored to the specific
calculation to be performed. Therefore, at this point in time, OpenMC is not
distributed with any pre-existing multi-group cross section libraries.
However, if the user has obtained or generated their own library, the user
should set the :envvar:`OPENMC_MG_CROSS_SECTIONS` environment variable
to the absolute path of the file library expected to used most frequently.
.. _NJOY: http://t2.lanl.gov/nis/codes.shtml
.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html
.. _NEA: http://www.oecd-nea.org

302
docs/source/usersguide/mgxs_library.rst Normal file
View file

@ -0,0 +1,302 @@
.. _usersguide_mgxs_library:
========================================
Multi-Group Cross Section Library Format
========================================
OpenMC can be run in continuous-energy mode or multi-group mode, provided the
nuclear data is available. In continuous-energy mode, the
``cross_sections.xml`` file contains necessary meta-data for each data set,
including the name and a file system location where the complete library
can be found. In multi-group mode, this ``cross_sections.xml`` file contains
this same meta-data describing the nuclide or material, but also contains the
group-wise nuclear data. This portion of the manual describes the format of
the multi-group data library required to be used in the ``cross_sections.xml``
file.
Similar to the other input file types, the multi-group library is provided in
the XML_ format. This library must provide some meta-data about the library
itself (such as the number of groups and the group structure, etc.) as well as
the actual cross section data itself for each of the necessary nuclides or
materials.
.. _XML: http://www.w3.org/XML/
------------------------------------------------
MGXS Library Specification -- cross_sections.xml
------------------------------------------------
The multi-group library meta-data is contained within the groups_,
group_structure_, and inverse_velocities_ elements.
The actual multi-group data itself is contained within the xsdata_ element.
.. _groups:
``<groups>`` Element
----------------------------------
The ``<groups>`` element has no attributes and simply provides the number of
energy groups contained within the library.
*Default*: None, this must be provided.
.. _group_structure:
``<group_structure>`` Element
-----------------------------
The ``<group_structure>`` element has no attributes and should be provided as a
monotonically increasing list of bounding energies, in MeV, for a number of
groups. To provide proper energy boundaries, the length of the data within the
``<group_structure>`` element should be one more than the number of groups in
the problem. For example, a two-group problem could be specified as:
.. code-block:: xml
<group_structure> 0.0 0.625E-6 20.0 </group_structure>
*Default*: None, this must be provided.
.. _inverse_velocities:
``<inverse_velocities>`` Element
--------------------------------
The ``<inverse_velocities>`` element optionally indicates the average
inverse velocity corresponding to each of the groups in the problem.
This element should therefore be an array with a length which matches the
number of groups set in the groups_ element.
*Default*: Should this be needed by the presence of an ``inverse-velocity``
score in the ``tallies.xml`` file and not provided in this element, OpenMC
will simply convert the group mid-point energy to an inverse of the velocity
and use this information for tallying.
.. _xsdata:
``<xsdata>`` Element
--------------------
The ``<xsdata>`` element contains the nuclide or material-specific meta-data as
well as the actual cross section data. The following are the
attributes/sub-elements required to describe the meta-data:
:name:
The name of the microscopic or macroscopic data set. An extension to the
name must be provided (e.g., the ``.300K`` in ``UO2.300K``). The name and
extension together must be twelve or less characters in length. This
extension must follow a period and be five characters or less in length.
similar to the equivalent in the continuous-energy ``cross_sections.xml``
file, is used to denote variants of the particular nuclide or material of
interest (i.e. the ``UO2`` data in this example could have been generated
at a temperature of 300K).
*Default*: None, this must be provided.
:alias:
An alternative name to use for the microscopic or macroscopic data set.
*Default*: If no alias is provided, it will adopt the value of ``name``.
:kT:
The temperature times Boltzmann's constant (in units of MeV) at which the
data was generated.
*Default*: Room temperature, 2.53E-8 MeV
:fissionable:
This element states whether or not the data in question is fissionable.
Accepted values are "true" or "false".
*Default*: None, this element must be provided.
:representation:
This element provides the method used to generate and represent the
multi-group cross sections. That is, whether they were generated with
scalar flux weighting (or reduced to an equivalent representation)
and thus are angle-independent, or if the data was generated with angular
dependent fluxes and thus the data is angle-dependent. The options are
either "isotropic" or "angle".
*Default*: "isotropic"
:num_azimuthal:
This element provides the number of equal width angular bins that the
azimuthal angular domain is subdivided in the case of angle-dependent
cross sections (i.e., "angle" is passed to the ``representation`` element).
Note that these bins are equal in azimuthal angle widths, not equal in the
cosine of the azimuthal angle widths.
*Default*: If ``representation`` is "angle", this must be provided. This
parameter is not used for other ``representation`` types.
:num_polar:
This element provides the number of equal width angular bins that the
polar angular domain is subdivided in the case of angle-dependent
cross sections (i.e., "angle" is passed to the ``representation`` element).
Note that these bins are equal in polar angle widths, not equal in the
cosine of the polar angle widths.
*Default*: If ``representation`` is "angle", this must be provided. This
parameter is not used for other ``representation`` types.
:scatt_type:
This element provides the representation of the angular distribution
associated with each group-to-group transfer probability. The options are
either "legendre", "histogram", or "tabular".
The "legendre" option means the angular distribution has been
expanded via Legendre polynomials of the order provided in the "order"
element.
The "histogram" option means the angular distribution is provided in
an equi-width histogram format with a number of bins as provided in the
"order" element. This is useful when the angular distribution was
obtained from a Monte Carlo tally and thus is natively in the histogram
format.
The "tabular" option means the angular distribution is provided in an
equi-spaced point-wise representation.
*Default*: "legendre"
:order:
This element provides either the Legendre order, number of bins, or number
of points used to describe the angular distribution associated with each
group-to-group transfer probability. The specific meaning of this bin
depends upon the value of ``scatt_type`` as discussed above.
*Default*: None, this element must be provided.
:tabular_legendre:
This optional element is used to set how the Legendre scattering kernel, if
provided via the ``scatt_type`` element above, is represented and thus used
during the scattering process. Specifically, the options are to either
convert the Legendre expansion to a tabular representation or leave it as
a set of Legendre coefficients. Converting to a tabular representation will
cost memory but is likely to decrease runtime compared to leaving as a
set of Legendre coefficients. This element has the following
attributes/sub-elements:
:enable:
This attribute/sub-element denotes whether or not the conversion to the
tabular format should be performed or not. A value of "true" means
the conversion should be performed, "false" means it should not.
*Default*: "true"
:num_points:
If the conversion is to take place the number of tabular points is
required. This attribute/sub-element allows the user to set the desired
number of points.
*Default*: 33
The following attributes/sub-elements are the cross section values to
be used during the transport process.
:total:
This element requires the group-wise total cross section ordered by
increasing group index (i.e., fast to thermal). If ``representation`` is
"isotropic", then the length of this list should equal the number of
groups described in the ``groups`` element. If ``representation`` is
"angle", then the length of this list should equal the number of groups
times the number of azimuthal angles times the number of polar angles,
with the inner-dimension being groups, intermediate-dimension being
azimuthal angles and outer-dimension being the polar angles.
*Default*: If not provided, it will be determined by summing the
absorption and scattering cross sections.
:absorption:
This element requires the group-wise absorption cross section ordered by
increasing group index (i.e., fast to thermal). If ``representation`` is
"isotropic", then the length of this list should equal the number of
groups described in the ``groups`` element. If ``representation`` is
"angle", then the length of this list should equal the number of groups
times the number of azimuthal angles times the number of polar angles,
with the inner-dimension being groups, intermediate-dimension being
azimuthal angles and outer-dimension being the polar angles.
*Default*: None, this must be provided.
:scatter:
This element requires the scattering moment matrices presented with the
columns representing incoming group and rows representing the outgoing
group. That is, down-scatter will be above the diagonal of the resultant
matrix. This matrix is repeated for every Legendre order (in order of
increasing orders) if ``scatt_type`` is "legendre"; otherwise, this
matrix is repeated for every bin of the histogram or tabular
representation. Finally, if ``representation`` is "angle", the above
is repeated for every azimuthal angle and every polar angle, in that
order.
*Default*: None, this must be provided.
:multiplicity:
This element provides the ratio of neutrons produced in scattering
collisions to the neutrons which undergo scattering collisions; that is,
the multiplicity provides the code with a scaling factor to account for
neutrons being produced in (n,xn) reactions. This information is assumed
isotropic and therefore does not need to be repeated for every Legendre
moment or histogram/tabular bin. This matrix follows the same arrangement
as described for the ``scatter`` element, with the exception of the
data needed to provide the scattering type information.
*Default*: Multiplicities of 1.0 are assumed (i.e., (n,xn) reactions are
neglected).
The following fission-specific data are only needed should ``fissionable``
be "true".
:fission:
This element requires the group-wise fission cross section ordered by
increasing group index (i.e., fast to thermal). If ``representation`` is
"isotropic", then the length of this list should equal the number of
groups described in the ``groups`` element. If ``representation`` is
"angle", then the length of this list should equal the number of groups
times the number of azimuthal angles times the number of polar angles,
with the inner-dimension being groups, intermediate-dimension being
azimuthal angles and outer-dimension being the polar angles.
*Default*: None, this is required only if fission tallies are
requested and the material is fissionable.
:kappa_fission:
This element requires the group-wise kappa-fission cross section ordered by
increasing group index (i.e., fast to thermal). If ``representation`` is
"isotropic", then the length of this list should equal the number of
groups described in the ``groups`` element. If ``representation`` is
"angle", then the length of this list should equal the number of groups
times the number of azimuthal angles times the number of polar angles,
with the inner-dimension being groups, intermediate-dimension being
azimuthal angles and outer-dimension being the polar angles.
*Default*: None, this is required only if kappa_fission tallies are
requested and the material is fissionable.
:chi:
This element requires the group-wise fission spectra ordered by
increasing group index (i.e., fast to thermal). This element should be
used if making the common approximation that the fission spectra does
not depend on incoming energy. If the user does not wish to make this
approximation, then this should not be provided and this information
included in the ``nu_fission`` element instead. If ``representation`` is
"isotropic", then the length of this list should equal the number of
groups described in the ``groups`` element. If ``representation`` is
"angle", then the length of this list should equal the number of groups
times the number of azimuthal angles times the number of polar angles,
with the inner-dimension being groups, intermediate-dimension being
azimuthal angles and outer-dimension being the polar angles.
*Default*: None, either this element is provided or ``nu_fission`` is
provided in fission matrix form, or the material is not fissionable.
:nu_fission:
This element provides either the group-wise fission production cross
section vector (i.e., if ``chi`` is provided), or is the group-wise fission
production matrix. If providing the vector, it should be ordered the same
as the ``fission`` data. If providing the matrix, it should be ordered
the same as the ``multiplicity`` matrix.
*Default*: None, either this element must be provided if the material
is fissionable.

3
docs/source/usersguide/output/particle_restart.rst
View file

@ -46,7 +46,8 @@ The current revision of the particle restart file format is 1.
**/energy** (*double*)
Energy of the particle in MeV.
Energy of the particle in MeV for continuous-energy mode, or the energy
group of the particle for multi-group mode.
**/xyz** (*double[3]*)

5
docs/source/usersguide/output/source.rst
View file

@ -15,5 +15,6 @@ is that documented here.
**/source_bank** (Compound type)
Source bank information for each particle. The compound type has fields
``wgt``, ``xyz``, ``uvw``, and ``E`` which represent the weight, position,
direction, and energy of the source particle, respectively.
``wgt``, ``xyz``, ``uvw``, ``E``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.

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

@ -4,7 +4,7 @@
State Point File Format
=======================
The current revision of the statepoint file format is 14.
The current revision of the statepoint file format is 15.
**/filetype** (*char[]*)
@ -39,6 +39,12 @@ The current revision of the statepoint file format is 14.
Pseudo-random number generator seed.
**/run_CE** (*int*)
Flag to denote continuous-energy or multi-group mode. A value of 1
indicates a continuous-energy run while a value of 0 indicates a
multi-group run.
**/run_mode** (*char[]*)
Run mode used. A value of 1 indicates a fixed-source run and a value of 2
@ -251,5 +257,6 @@ if (run_mode == 'k-eigenvalue' and source_present > 0)
**/source_bank** (Compound type)
Source bank information for each particle. The compound type has fields
``wgt``, ``xyz``, ``uvw``, and ``E`` which represent the weight,
position, direction, and energy of the source particle, respectively.
``wgt``, ``xyz``, ``uvw``, ``E``, ``g``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.

180
examples/python/pincell_multigroup/build-xml.py Normal file
View file

@ -0,0 +1,180 @@
import openmc
import openmc.mgxs
from openmc.source import Source
from openmc.stats import Box
import numpy as np
###############################################################################
# Simulation Input File Parameters
###############################################################################
# OpenMC simulation parameters
batches = 100
inactive = 10
particles = 1000
###############################################################################
# Exporting to OpenMC mg_cross_sections.xml File
###############################################################################
# Instantiate the energy group data
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.order = 0
uo2_xsdata.total = np.array([0.1779492, 0.3298048, 0.4803882, 0.5543674,
0.3118013, 0.3951678, 0.5644058])
uo2_xsdata.absorption = np.array([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
3.0020E-02, 1.1126E-01, 2.8278E-01])
scatter = [[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.4525650, 0.0055664, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]]
uo2_xsdata.scatter = np.array(scatter[:][:])
uo2_xsdata.fission = np.array([7.21206E-03, 8.19301E-04, 6.45320E-03,
1.85648E-02, 1.78084E-02, 8.30348E-02,
2.16004E-01])
uo2_xsdata.nu_fission = np.array([2.005998E-02, 2.027303E-03, 1.570599E-02,
4.518301E-02, 4.334208E-02, 2.020901E-01,
5.257105E-01])
uo2_xsdata.chi = np.array([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.order = 0
h2o_xsdata.total = np.array([0.15920605, 0.412969593, 0.59030986, 0.58435,
0.718, 1.2544497, 2.650379])
h2o_xsdata.absorption = np.array([6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02])
scatter = [[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
[0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
[0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
[0.0000000, 0.0000000, 0.0000000, 0.0910284, 0.4155100, 0.0637320, 0.0121390],
[0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]]
h2o_xsdata.scatter = np.array(scatter)
mg_cross_sections_file = openmc.MGXSLibraryFile(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata,h2o_xsdata])
mg_cross_sections_file.export_to_xml()
###############################################################################
# Exporting to OpenMC materials.xml File
###############################################################################
# Instantiate some Macroscopic Data
uo2_data = openmc.Macroscopic('UO2', '300K')
h2o_data = openmc.Macroscopic('LWTR', '300K')
# Instantiate some Materials and register the appropriate Macroscopic objects
uo2 = openmc.Material(material_id=1, name='UO2 fuel')
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
water = openmc.Material(material_id=2, name='Water')
water.set_density('macro', 1.0)
water.add_macroscopic(h2o_data)
# Instantiate a MaterialsFile, register all Materials, and export to XML
materials_file = openmc.MaterialsFile()
materials_file.default_xs = '300K'
materials_file.add_materials([uo2, water])
materials_file.export_to_xml()
###############################################################################
# Exporting to OpenMC geometry.xml File
###############################################################################
# Instantiate ZCylinder surfaces
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.54, name='Fuel OR')
left = openmc.XPlane(surface_id=4, x0=-0.63, name='left')
right = openmc.XPlane(surface_id=5, x0=0.63, name='right')
bottom = openmc.YPlane(surface_id=6, y0=-0.63, name='bottom')
top = openmc.YPlane(surface_id=7, y0=0.63, name='top')
left.boundary_type = 'reflective'
right.boundary_type = 'reflective'
top.boundary_type = 'reflective'
bottom.boundary_type = 'reflective'
# Instantiate Cells
fuel = openmc.Cell(cell_id=1, name='cell 1')
moderator = openmc.Cell(cell_id=2, name='cell 2')
# Use surface half-spaces to define regions
fuel.region = -fuel_or
moderator.region = +fuel_or & +left & -right & +bottom & -top
# Register Materials with Cells
fuel.fill = uo2
moderator.fill = water
# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
# Register Cells with Universe
root.add_cells([fuel, moderator])
# Instantiate a Geometry and register the root Universe
geometry = openmc.Geometry()
geometry.root_universe = root
# Instantiate a GeometryFile, register Geometry, and export to XML
geometry_file = openmc.GeometryFile()
geometry_file.geometry = geometry
geometry_file.export_to_xml()
###############################################################################
# Exporting to OpenMC settings.xml File
###############################################################################
# Instantiate a SettingsFile, set all runtime parameters, and export to XML
settings_file = openmc.SettingsFile()
settings_file.energy_mode = "multi-group"
settings_file.cross_sections = "./mg_cross_sections.xml"
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.source = Source(space=Box([-0.63, -0.63, -1.], [0.63, 0.63, 1.]))
###############################################################################
# Exporting to OpenMC tallies.xml File
###############################################################################
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [100, 100, 1]
mesh.lower_left = [-0.63, -0.63, -1.e50]
mesh.upper_right = [0.63, 0.63, 1.e50]
# Instantiate some tally Filters
energy_filter = openmc.Filter(type='energy',
bins=[1E-11, 0.0635E-6, 10.0E-6, 1.0E-4, 1.0E-3,
0.5, 1.0, 20.0])
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
# Instantiate the Tally
tally = openmc.Tally(tally_id=1, name='tally 1')
tally.add_filter(energy_filter)
tally.add_filter(mesh_filter)
tally.add_score('flux')
tally.add_score('fission')
tally.add_score('nu-fission')
# Instantiate a TalliesFile, register all Tallies, and export to XML
tallies_file = openmc.TalliesFile()
tallies_file.add_mesh(mesh)
tallies_file.add_tally(tally)
tallies_file.export_to_xml()

16
examples/xml/pincell_multigroup/geometry.xml Normal file
View file

@ -0,0 +1,16 @@
<geometry>
<surface coeffs="0. 0. 0.540" id="1" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="20" type="x-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="21" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="22" type="y-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="23" type="y-plane" />
<cell id="1" material="1" region=" -1" />
<cell id="2" material="7" region="1 20 -21 22 -23" />
</geometry>

54
examples/xml/pincell_multigroup/materials.xml Normal file
View file

@ -0,0 +1,54 @@
<?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" />
<macroscopic name="UO2"/>
</material>
<!-- 4.3% MOX -->
<material id="2">
<density units="macro" value="1.0" />
<macroscopic name="MOX1"/>
</material>
<!-- 7.0 MOX -->
<material id="3">
<density units="macro" value="1.0" />
<macroscopic name="MOX2"/>
</material>
<!-- 8.0% MOX -->
<material id="4">
<density units="macro" value="1.0" />
<macroscopic name="MOX3"/>
</material>
<!-- Fission Chamber -->
<material id="5">
<density units="macro" value="1.0" />
<macroscopic name="FC"/>
</material>
<!-- Guide Tube -->
<material id="6">
<density units="macro" value="1.0" />
<macroscopic name="GT"/>
</material>
<!-- Water -->
<material id="7">
<density units="macro" value="1.0" />
<macroscopic name="LWTR"/>
</material>
<!-- Control Rod -->
<material id="8">
<density units="macro" value="1.0" />
<macroscopic name="CR"/>
</material>
</materials>

383
examples/xml/pincell_multigroup/mg_cross_sections.xml Normal file
View file

@ -0,0 +1,383 @@
<?xml version="1.0"?>
<library>
<!-- Before getting to the data, set common information -->
<groups> 7 </groups>
<group_structure>
1E-11 0.0635E-6 10.0E-6 1.0E-4 1.0E-3 0.5 1.0 20.0
</group_structure>
<!--
Move on to the data. Each <xsdata> has a unique id and label.
-->
<xsdata>
<!-- Meta data for this data -->
<name>UO2.300K</name>
<alias>UO2.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
<!-- Optional (default is isotropic) -->
<representation>isotropic</representation>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
8.0248E-03 3.7174E-03 2.6769E-02 9.6236E-02 3.0020E-02 1.1126E-01 2.8278E-01
</absorption>
<nu_fission>
2.005998E-02 2.027303E-03 1.570599E-02 4.518301E-02 4.334208E-02 2.020901E-01 5.257105E-01
</nu_fission>
<chi>
5.8791E-01 4.1176E-01 3.3906E-04 1.1761E-07 0.0000E+00 0.0000E+00 0.0000E+00
</chi>
<fission>
7.21206E-03 8.19301E-04 6.45320E-03 1.85648E-02 1.78084E-02 8.30348E-02 2.16004E-01
</fission>
<!-- units of MeV/cm -->
<!-- If no kappa fission tallies, this is not needed; it will not be loaded
if there is no kappa fission scores anyways -->
<k_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</k_fission>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
0.1275370 0.0423780 0.0000094 0.0000000 0.0000000 0.0000000 0.0000000
0.0000000 0.3244560 0.0016314 0.0000000 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.4509400 0.0026792 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000 0.4525650 0.0055664 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000 0.0001253 0.2714010 0.0102550 0.0000000
0.0000000 0.0000000 0.0000000 0.0000000 0.0012968 0.2658020 0.0168090
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0085458 0.2730800
</scatter>
<!-- If total is not provided, it will be calculated.
However, in the C5G7 problems, we want to use a transport-corrected value
so we dont want to have it be calculated -->
<total>
0.1779492 0.3298048 0.4803882 0.5543674000000001 0.3118013 0.39516779999999996 0.5644058
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>MOX1.300K</name>
<alias>MOX1.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
8.4339E-03 3.7577E-03 2.7970E-02 1.0421E-01 1.3994E-01 4.0918E-01 4.0935E-01
</absorption>
<!--
Since chi_vector is false, this will be a matrix
Matrix is g_in, g_out.
This is to show that you can either use a chi vector + nu_fission vector,
like in the UO2 data, or a nu_fission matrix like here.
-->
<nu_fission>
1.27888062E-02 8.95701528E-03 7.37557218E-06 2.55837033E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
1.49041240E-03 1.04385401E-03 8.59552023E-07 2.98153464E-10 0.00000000E+00 0.00000000E+00 0.00000000E+00
9.56411400E-03 6.69850756E-03 5.51582469E-06 1.91327830E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
3.84928781E-02 2.69596154E-02 2.21996483E-05 7.70040890E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
1.80629998E-02 1.26509513E-02 1.04173100E-05 3.61346022E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
3.91930789E-01 2.74500216E-01 2.26034688E-04 7.84048241E-08 0.00000000E+00 0.00000000E+00 0.00000000E+00
4.19762096E-01 2.93992687E-01 2.42085585E-04 8.39724109E-08 0.00000000E+00 0.00000000E+00 0.00000000E+00
</nu_fission>
<fission>
7.62704E-03 8.76898E-04 5.69835E-03 2.28872E-02 1.07635E-02 2.32757E-01 2.48968E-01
</fission>
<!-- units of MeV/cm -->
<k_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</k_fission>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
1.27537000E-01 4.23780000E-02 9.43740000E-06 5.51630000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 3.24456000E-01 1.63140000E-03 3.14270000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 4.50940000E-01 2.67920000E-03 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 4.52565000E-01 5.56640000E-03 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 1.25250000E-04 2.71401000E-01 1.02550000E-02 1.00210000E-08
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 1.29680000E-03 2.65802000E-01 1.68090000E-02
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 8.54580000E-03 2.73080000E-01
</scatter>
<total>
0.1783583429163 0.3298451031427 0.4815892 0.5623414 0.421721260021 0.6930878 0.6909757999999999
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>MOX2.300K</name>
<alias>MOX2.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
0.0090657 0.0042967 0.032881 0.12203 0.18298 0.56846 0.58521
</absorption>
<!--
Since chi_vector is false, this will be a matrix
Matrix is g_in, g_out !!! Need to get these values looking right to match
output of a nu-fission tally with <energy> filter above <energyout> filter
-->
<nu_fission>
1.40004593E-02 9.80563205E-03 8.07435789E-06 2.80075866E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
2.26856185E-03 1.58885378E-03 1.30832709E-06 4.53820413E-10 0.00000000E+00 0.00000000E+00 0.00000000E+00
1.41886199E-02 9.93741584E-03 8.18287404E-06 2.83839974E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
5.54788444E-02 3.88562347E-02 3.19958106E-05 1.10984111E-08 0.00000000E+00 0.00000000E+00 0.00000000E+00
2.69085702E-02 1.88462058E-02 1.55187355E-05 5.38299559E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
5.45687127E-01 3.82187973E-01 3.14709185E-04 1.09163414E-07 0.00000000E+00 0.00000000E+00 0.00000000E+00
6.13307712E-01 4.29548032E-01 3.53707392E-04 1.22690752E-07 0.00000000E+00 0.00000000E+00 0.00000000E+00
</nu_fission>
<fission>
0.00825446 0.00132565 0.00842156 0.032873 0.0159636 0.323794 0.362803
</fission>
<!-- units of MeV/cm -->
<k_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</k_fission>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
1.30457000E-01 4.17920000E-02 8.51050000E-06 5.13290000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 3.28428000E-01 1.64360000E-03 2.20170000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 4.58371000E-01 2.53310000E-03 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 4.63709000E-01 5.47660000E-03 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 1.76190000E-04 2.82313000E-01 8.72890000E-03 9.00160000E-09
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2.27600000E-03 2.49751000E-01 1.31140000E-02
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 8.86450000E-03 2.59529000E-01
</scatter>
<total>
0.1813232156329 0.3343683022017 0.4937851 0.5912156 0.47419809900160004 0.833601 0.8536035
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>MOX3.300K</name>
<alias>MOX3.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
9.48620000E-03 4.65560000E-03 3.62400000E-02 1.32720000E-01 2.08400000E-01 6.58700000E-01 6.90170000E-01
</absorption>
<!--
Since chi_vector is false, this will be a matrix
Matrix is g_in, g_out !!! Need to get these values looking right to match
output of a nu-fission tally with <energy> filter above <energyout> filter
-->
<nu_fission>
1.48071013E-02 1.03705874E-02 8.53956516E-06 2.96212546E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
2.78640474E-03 1.95154023E-03 1.60697792E-06 5.57413653E-10 0.00000000E+00 0.00000000E+00 0.00000000E+00
1.73304404E-02 1.21378819E-02 9.99482763E-06 3.46691346E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
6.59928975E-02 4.62200600E-02 3.80594850E-05 1.32017225E-08 0.00000000E+00 0.00000000E+00 0.00000000E+00
3.25131926E-02 2.27715674E-02 1.87510386E-05 6.50418701E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
6.32002662E-01 4.42641588E-01 3.64489161E-04 1.26430632E-07 0.00000000E+00 0.00000000E+00 0.00000000E+00
7.28595687E-01 5.10293344E-01 4.20196380E-04 1.45753838E-07 0.00000000E+00 0.00000000E+00 0.00000000E+00
</nu_fission>
<fission>
8.67209000E-03 1.62426000E-03 1.02716000E-02 3.90447000E-02 1.92576000E-02 3.74888000E-01 4.30599000E-01
</fission>
<!-- units of MeV/cm -->
<k_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</k_fission>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
1.31504000E-01 4.20460000E-02 8.69720000E-06 5.19380000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 3.30403000E-01 1.64630000E-03 2.60060000E-09 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 4.61792000E-01 2.47490000E-03 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 4.68021000E-01 5.43300000E-03 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 1.85970000E-04 2.85771000E-01 8.39730000E-03 8.92800000E-09
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2.39160000E-03 2.47614000E-01 1.23220000E-02
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 8.96810000E-03 2.56093000E-01
</scatter>
<total>
1.83044902E-01 3.36704903E-01 5.00506900E-01 6.06174000E-01 5.02754279E-01 9.21027600E-01 9.55231100E-01
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>FC.300K</name>
<alias>FC.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>true</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
5.1132E-04 7.5813E-05 3.1643E-04 1.1675E-03 3.3977E-03 9.1886E-03 2.3244E-02
</absorption>
<nu_fission>
1.323401E-08 1.434500E-08 1.128599E-06 1.276299E-05 3.538502E-07 1.740099E-06 5.063302E-06
</nu_fission>
<chi>
5.8791E-01 4.1176E-01 3.3906E-04 1.1761E-07 0.0000E+00 0.0000E+00 0.0000E+00
</chi>
<fission>
4.79002E-09 5.82564E-09 4.63719E-07 5.24406E-06 1.45390E-07 7.14972E-07 2.08041E-06
</fission>
<!-- units of MeV/cm -->
<k_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</k_fission>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
6.61659000E-02 5.90700000E-02 2.83340000E-04 1.46220000E-06 2.06420000E-08 0.00000000E00 0.00000000E00
0.00000000E00 2.40377000E-01 5.24350000E-02 2.49900000E-04 1.92390000E-05 2.98750000E-06 4.21400000E-07
0.00000000E00 0.00000000E00 1.83425000E-01 9.22880000E-02 6.93650000E-03 1.07900000E-03 2.05430000E-04
0.00000000E00 0.00000000E00 0.00000000E00 7.90769000E-02 1.69990000E-01 2.58600000E-02 4.92560000E-03
0.00000000E00 0.00000000E00 0.00000000E00 3.73400000E-05 9.97570000E-02 2.06790000E-01 2.44780000E-02
0.00000000E00 0.00000000E00 0.00000000E00 0.00000000E00 9.17420000E-04 3.16774000E-01 2.38760000E-01
0.00000000E00 0.00000000E00 0.00000000E00 0.00000000E00 0.00000000E00 4.97930000E-02 1.09910000E00
</scatter>
<total>
1.26032048E-01 2.93160367E-01 2.84250824E-01 2.81025244E-01 3.34460185E-01 5.65640735E-01 1.17213908E00
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>GT.300K</name>
<alias>GT.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
5.11320000E-04 7.58010000E-05 3.15720000E-04 1.15820000E-03 3.39750000E-03 9.18780000E-03 2.32420000E-02
</absorption>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
6.61659000E-02 5.90700000E-02 2.83340000E-04 1.46220000E-06 2.06420000E-08 0.00000000E+00 0.00000000E+00
0.00000000E+00 2.40377000E-01 5.24350000E-02 2.49900000E-04 1.92390000E-05 2.98750000E-06 4.21400000E-07
0.00000000E+00 0.00000000E+00 1.83297000E-01 9.23970000E-02 6.94460000E-03 1.08030000E-03 2.05670000E-04
0.00000000E+00 0.00000000E+00 0.00000000E+00 7.88511000E-02 1.70140000E-01 2.58810000E-02 4.92970000E-03
0.00000000E+00 0.00000000E+00 0.00000000E+00 3.73330000E-05 9.97372000E-02 2.06790000E-01 2.44780000E-02
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 9.17260000E-04 3.16765000E-01 2.38770000E-01
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 4.97920000E-02 1.09912000E+00
</scatter>
<total>
1.26032043E-01 2.93160349E-01 2.84240290E-01 2.80960000E-01 3.34440033E-01 5.65640060E-01 1.17215400E+00
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>LWTR.300K</name>
<alias>LWTR.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
6.0105E-04 1.5793E-05 3.3716E-04 1.9406E-03 5.7416E-03 1.5001E-02 3.7239E-02
</absorption>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
0.0444777 0.1134000 0.0007235 0.0000037 0.0000001 0.0000000 0.0000000
0.0000000 0.2823340 0.1299400 0.0006234 0.0000480 0.0000074 0.0000010
0.0000000 0.0000000 0.3452560 0.2245700 0.0169990 0.0026443 0.0005034
0.0000000 0.0000000 0.0000000 0.0910284 0.4155100 0.0637320 0.0121390
0.0000000 0.0000000 0.0000000 0.0000714 0.1391380 0.5118200 0.0612290
0.0000000 0.0000000 0.0000000 0.0000000 0.0022157 0.6999130 0.5373200
0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1324400 2.4807000
</scatter>
<total>
0.15920605 0.41296959299999997 0.59030986 0.5843499999999999 0.7180000000000001 1.2544497000000001 2.650379
</total>
</xsdata>
<xsdata>
<!-- Meta data for this data -->
<name>CR.300K</name>
<alias>CR.300K</alias>
<kT> 2.53E-8 </kT> <!-- in MeV -->
<order>0</order>
<fissionable>false</fissionable>
<!-- The data itself, like tallies,
goes from low energies (groups) to high energies
-->
<absorption>
1.70490000E-03 8.36224000E-03 8.37901000E-02 3.97797000E-01 6.98763000E-01 9.29508000E-01 1.17836000E+00
</absorption>
<!-- for consistency must include nu-scatter -->
<!-- will be a matrix of (order+1) x g_in x g_out -->
<scatter>
1.70563000E-01 4.44012000E-02 9.83670000E-05 1.27786000E-07 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 4.71050000E-01 6.85480000E-04 3.91395000E-10 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 8.01859000E-01 7.20132000E-04 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 5.70752000E-01 1.46015000E-03 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 6.55562000E-05 2.07838000E-01 3.81486000E-03 3.69760000E-09
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 1.02427000E-03 2.02465000E-01 4.75290000E-03
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 3.53043000E-03 6.58597000E-01
</scatter>
<total>
2.16767595E-01 4.80097720E-01 8.86369232E-01 9.70009150E-01 9.10481420E-01 1.13775017E+00 1.84048743E+00
</total>
</xsdata>
</library>

28
examples/xml/pincell_multigroup/plots.xml Normal file
View file

@ -0,0 +1,28 @@
<?xml version="1.0"?>
<plots>
<plot>
<id>1</id>
<filename>mat</filename>
<color>material</color>
<origin>0 0 0</origin>
<width>1.26 1.26</width>
<type>slice</type>
<pixels>1000 1000 </pixels>
<col_spec id="1" rgb="255 0 0" />
<col_spec id="2" rgb="0 0 0" />
<col_spec id="3" rgb="0 255 0" />
<col_spec id="4" rgb="0 0 255" />
</plot>
<plot>
<id>2</id>
<filename>cell</filename>
<color>cell</color>
<origin>0 0 0</origin>
<width>1.26 1.26</width>
<type>slice</type>
<pixels>1000 1000 </pixels>
</plot>
</plots>

40
examples/xml/pincell_multigroup/settings.xml Normal file
View file

@ -0,0 +1,40 @@
<?xml version="1.0"?>
<settings>
<energy_mode>multi-group</energy_mode>
<!--
Define how many particles to run and for how many batches
in an eigenvalue calculation mode
-->
<eigenvalue>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!--
Start with uniformally distributed neutron source
with the default energy spectrum of a Maxwellian
and isotropic distribution.
-->
<source>
<space type="box">
<parameters>
-0.63 -0.63 -1E50
0.63 0.63 1E50
</parameters>
</space>
</source>
<output>
<cross_sections>true</cross_sections>
<summary>true</summary>
<tallies>true</tallies>
</output>
<survival_biasing>false</survival_biasing>
<cross_sections>./mg_cross_sections.xml</cross_sections>
</settings>

13
examples/xml/pincell_multigroup/tallies.xml Normal file
View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1" type="regular">
<dimension>100 100 1</dimension>
<lower_left>-0.63 -0.63 -1e+50</lower_left>
<upper_right>0.63 0.63 1e+50</upper_right>
</mesh>
<tally id="1" name="tally 1">
<filter bins="1e-11 6.35e-08 1e-05 0.0001 0.001 0.5 1.0 20.0" type="energy" />
<filter bins="1" type="mesh" />
<scores>flux fission nu-fission</scores>
</tally>
</tallies>

2
openmc/__init__.py
View file

@ -1,11 +1,13 @@
from openmc.element import *
from openmc.geometry import *
from openmc.nuclide import *
from openmc.macroscopic import *
from openmc.material import *
from openmc.plots import *
from openmc.settings import *
from openmc.surface import *
from openmc.universe import *
from openmc.mgxs_library import *
from openmc.mesh import *
from openmc.filter import *
from openmc.trigger import *

2
openmc/checkvalue.py
View file

@ -127,7 +127,7 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1):
# But first, have we exceeded the max depth?
if len(tree) > max_depth:
msg = 'Error setting {0}: Found an iterable at {1}, items '\
'in that iterable excceed the maximum depth of {2}' \
'in that iterable exceed the maximum depth of {2}' \
.format(name, ind_str, max_depth)
raise ValueError(msg)

85
openmc/macroscopic.py Normal file
View file

@ -0,0 +1,85 @@
from numbers import Integral
import sys
from openmc.checkvalue import check_type
if sys.version_info[0] >= 3:
basestring = str
class Macroscopic(object):
"""A Macroscopic object that can be used in a material.
Parameters
----------
name : str
Name of the macroscopic data, e.g. UO2
xs : str
Cross section identifier, e.g. 71c
Attributes
----------
name : str
Name of the nuclide, e.g. UO2
xs : str
Cross section identifier, e.g. 71c
"""
def __init__(self, name='', xs=None):
# Initialize class attributes
self._name = ''
self._xs = None
# Set the Material class attributes
self.name = name
if xs is not None:
self.xs = xs
def __eq__(self, other):
if isinstance(other, Macroscopic):
if self._name != other._name:
return False
elif self._xs != other._xs:
return False
else:
return True
elif isinstance(other, basestring) and other == self.name:
return True
else:
return False
def __ne__(self, other):
return not self == other
def __hash__(self):
return hash((self._name, self._xs))
def __repr__(self):
string = 'Nuclide - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
return string
@property
def name(self):
return self._name
@property
def xs(self):
return self._xs
@name.setter
def name(self, name):
check_type('name', name, basestring)
self._name = name
@xs.setter
def xs(self, xs):
check_type('cross-section identifier', xs, basestring)
self._xs = xs
def __repr__(self):
string = 'Macroscopic - {0}\n'.format(self._name)
string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self.xs)
return string

140
openmc/material.py
View file

@ -22,14 +22,15 @@ def reset_auto_material_id():
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum']
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum',
'macro']
# Constant for density when not needed
NO_DENSITY = 99999.
class Material(object):
"""A material composed of a collection of nuclides/elements that can be
"""A material composed of a collection of nuclides/elements that can be
assigned to a region of space.
Parameters
@ -49,7 +50,8 @@ class Material(object):
Density of the material (units defined separately)
density_units : str
Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3',
'atom/b-cm', 'atom/cm3', or 'sum'.
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
applies in the case of a multi-group calculation.
"""
@ -65,6 +67,10 @@ class Material(object):
# Values - tuple (nuclide, percent, percent type)
self._nuclides = OrderedDict()
# The single instance of Macroscopic data present in this material
# (only one is allowed, hence this is different than _nuclides, etc)
self._macroscopic = None
# An ordered dictionary of Elements (order affects OpenMC results)
# Keys - Element names
# Values - tuple (element, percent, percent type)
@ -128,6 +134,10 @@ class Material(object):
string += '{0: <16}'.format('\t{0}'.format(nuclide))
string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
if self._macroscopic is not None:
string += '{0: <16}\n'.format('\tMacroscopic Data')
string += '{0: <16}'.format('\t{0}'.format(self._macroscopic))
string += '{0: <16}\n'.format('\tElements')
for element in self._elements:
@ -149,6 +159,7 @@ class Material(object):
clone._density = self._density
clone._density_units = self._density_units
clone._nuclides = deepcopy(self._nuclides, memo)
clone._macroscopic = self._macroscopic
clone._elements = deepcopy(self._elements, memo)
clone._sab = deepcopy(self._sab, memo)
clone._convert_to_distrib_comps = self._convert_to_distrib_comps
@ -268,6 +279,11 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add a Nuclide to Material ID="{0}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(nuclide, (openmc.Nuclide, str)):
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-Nuclide value "{1}"'.format(self._id, nuclide)
@ -311,6 +327,64 @@ class Material(object):
if nuclide._name in self._nuclides:
del self._nuclides[nuclide._name]
def add_macroscopic(self, macroscopic):
"""Add a macroscopic to the material
Parameters
----------
macroscopic : str or Macroscopic
Macroscopic to add
"""
# Ensure no nuclides, elements, or sab are added since these would be
# incompatible with macroscopics
if self._nuclides or self._elements or self._sab:
msg = 'Unable to add a Macroscopic data set to Material ID="{0}" ' \
'with a macroscopic value "{1}" as an incompatible data ' \
'member (i.e., nuclide, element, or S(a,b) table) ' \
'has already been added'.format(self._id, macroscopic)
raise ValueError(msg)
if not isinstance(macroscopic, (openmc.Macroscopic, basestring)):
msg = 'Unable to add a Macroscopic to Material ID="{0}" with a ' \
'non-Macroscopic value "{1}"'.format(self._id, macroscopic)
raise ValueError(msg)
if isinstance(macroscopic, openmc.Macroscopic):
# Copy this Macroscopic to separate it from the Macroscopic in
# other Materials
macroscopic = deepcopy(macroscopic)
else:
macroscopic = openmc.Macroscopic(macroscopic)
if self._macroscopic is None:
self._macroscopic = macroscopic
else:
msg = 'Unable to add a Macroscopic to Material ID="{0}", ' \
'Only One Macroscopic allowed per ' \
'Material!'.format(self._id, macroscopic)
raise ValueError(msg)
def remove_macroscopic(self, macroscopic):
"""Remove a macroscopic from the material
Parameters
----------
macroscopic : Macroscopic
Macroscopic to remove
"""
if not isinstance(macroscopic, openmc.Macroscopic):
msg = 'Unable to remove a Macroscopic "{0}" in Material ID="{1}" ' \
'since it is not a Macroscopic'.format(self._id, macroscopic)
raise ValueError(msg)
# If the Material contains the Macroscopic, delete it
if macroscopic._name == self._macroscopic.name:
self._macroscopic = None
def add_element(self, element, percent, percent_type='ao'):
"""Add a natural element to the material
@ -325,6 +399,11 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add an Element to Material ID="{0}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(element, openmc.Element):
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-Element value "{1}"'.format(self._id, element)
@ -371,6 +450,11 @@ class Material(object):
"""
if self._macroscopic is not None:
msg = 'Unable to add an S(a,b) table to Material ID="{0}" as a ' \
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if not isinstance(name, basestring):
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string table name "{1}"'.format(self._id, name)
@ -427,6 +511,15 @@ class Material(object):
return xml_element
def _get_macroscopic_xml(self, macroscopic):
xml_element = ET.Element("macroscopic")
xml_element.set("name", macroscopic._name)
if macroscopic.xs is not None:
xml_element.set("xs", macroscopic.xs)
return xml_element
def _get_element_xml(self, element, distrib=False):
xml_element = ET.Element("element")
xml_element.set("name", str(element[0]._name))
@ -482,14 +575,19 @@ class Material(object):
subelement.set("units", self._density_units)
if not self._convert_to_distrib_comps:
# Create nuclide XML subelements
subelements = self._get_nuclides_xml(self._nuclides)
for subelement in subelements:
element.append(subelement)
if self._macroscopic is None:
# Create nuclide XML subelements
subelements = self._get_nuclides_xml(self._nuclides)
for subelement in subelements:
element.append(subelement)
# Create element XML subelements
subelements = self._get_elements_xml(self._elements)
for subelement in subelements:
# Create element XML subelements
subelements = self._get_elements_xml(self._elements)
for subelement in subelements:
element.append(subelement)
else:
# Create macroscopic XML subelements
subelement = self._get_macroscopic_xml(self._macroscopic)
element.append(subelement)
else:
@ -516,15 +614,21 @@ class Material(object):
subsubelement = ET.SubElement(subelement, "otf_file_path")
subsubelement.text = self._distrib_otf_file
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides, distrib=True)
for subelement_nuc in subelements:
subelement.append(subelement_nuc)
if self._macroscopic is None:
# Create nuclide XML subelements
subelements = self.get_nuclides_xml(self._nuclides, distrib=True)
for subelement_nuc in subelements:
subelement.append(subelement_nuc)
# Create element XML subelements
subelements = self._get_elements_xml(self._elements, distrib=True)
for subelement_ele in subelements:
subelement.append(subelement_ele)
# Create element XML subelements
subelements = self._get_elements_xml(self._elements, distrib=True)
for subsubelement in subelements:
subelement.append(subsubelement)
else:
# Create macroscopic XML subelements
subsubelement = self._get_macroscopic_xml(self._macroscopic,
distrib=True)
subelement.append(subsubelement)
if len(self._sab) > 0:
for sab in self._sab:

4
openmc/mgxs/groups.py
View file

@ -24,7 +24,7 @@ class EnergyGroups(object):
----------
group_edges : Iterable of Real
The energy group boundaries [MeV]
num_group : Integral
num_groups : Integral
The number of energy groups
"""
@ -54,7 +54,7 @@ class EnergyGroups(object):
def __eq__(self, other):
if not isinstance(other, EnergyGroups):
return False
elif self.group_edges != other.group_edges:
elif (self.group_edges != other.group_edges).all():
return False
else:
return True

722
openmc/mgxs_library.py Normal file
View file

@ -0,0 +1,722 @@
from collections import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import warnings
import sys
if sys.version_info[0] >= 3:
basestring = str
import numpy as np
import openmc
from openmc.mgxs import EnergyGroups
from openmc.checkvalue import check_type, check_value, check_greater_than, \
check_iterable_type
from openmc.clean_xml import *
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATIONS = ['isotropic', 'angle']
def ndarray_to_string(arr):
"""Converts a numpy ndarray in to a join with spaces between entries
similar to ' '.join(map(str,arr)) but applied to all sub-dimensions.
Parameters
----------
arr : ndarray
Array to combine in to a string
Returns
-------
text : str
String representation of array in arr
"""
shape = arr.shape
ndim = arr.ndim
tab = ' '
indent = '\n' + tab + tab
text = indent
if ndim == 1:
text += tab
for i in range(shape[0]):
text += '{:.7E} '.format(arr[i])
text += indent
elif ndim == 2:
for i in range(shape[0]):
text += tab
for j in range(shape[1]):
text += '{:.7E} '.format(arr[i, j])
text += indent
elif ndim == 3:
for i in range(shape[0]):
for j in range(shape[1]):
text += tab
for k in range(shape[2]):
text += '{:.7E} '.format(arr[i, j, k])
text += indent
elif ndim == 4:
for i in range(shape[0]):
for j in range(shape[1]):
for k in range(shape[2]):
text += tab
for l in range(shape[3]):
text += '{:.7E} '.format(arr[i, j, k, l])
text += indent
elif ndim == 5:
for i in range(shape[0]):
for j in range(shape[1]):
for k in range(shape[2]):
for l in range(shape[3]):
text += tab
for m in range(shape[4]):
text += '{:.7E} '.format(arr[i, j, k, l, m])
text += indent
return text
class XSdata(object):
"""A multi-group cross section data set providing all the
multi-group data necessary for a multi-group OpenMC calculation.
Parameters
----------
name : str, optional
Name of the mgxs data set.
representation : {'isotropic', 'angle'}
Method used in generating the MGXS (isotropic or angle-dependent flux
weighting). Defaults to 'isotropic'
Attributes
----------
name : str
Unique identifier for the xsdata object
alias : str
Separate unique identifier for the xsdata object
kT : float
Temperature (in units of MeV) of this data set.
energy_groups : EnergyGroups
Energy group structure
fissionable : boolean
Whether or not this is a fissionable data set.
scatt_type : {'legendre', 'histogram', or 'tabular'}
Angular distribution representation (legendre, histogram, or tabular)
order : int
Either the Legendre order, number of bins, or number of points used to
describe the angular distribution associated with each group-to-group
transfer probability.
tabular_legendre : dict
Set how to treat the Legendre scattering kernel (tabular or leave in
Legendre polynomial form). Dict contains two keys: 'enable' and
'num_points'. 'enable' is a boolean and 'num_points' is the
number of points to use, if 'enable' is True.
"""
def __init__(self, name, energy_groups, representation="isotropic"):
# Initialize class attributes
self._name = name
self._energy_groups = energy_groups
self._representation = representation
self._alias = None
self._kT = None
self._fissionable = False
self._scatt_type = 'legendre'
self._order = None
self._tabular_legendre = None
self._num_polar = None
self._num_azimuthal = None
self._total = None
self._absorption = None
self._scatter = None
self._multiplicity = None
self._fission = None
self._nu_fission = None
self._k_fission = None
self._chi = None
self._use_chi = None
@property
def name(self):
return self._name
@property
def energy_groups(self):
return self._energy_groups
@property
def representation(self):
return self._representation
@property
def alias(self):
return self._alias
@property
def kT(self):
return self._kT
@property
def scatt_type(self):
return self._scatt_type
@property
def order(self):
return self._order
@property
def tabular_legendre(self):
return self._tabular_legendre
@property
def num_polar(self):
return self._num_polar
@property
def num_azimuthal(self):
return self._num_azimuthal
@property
def total(self):
return self._total
@property
def absorption(self):
return self._absorption
@property
def scatter(self):
return self._scatter
@property
def multiplicity(self):
return self._multiplicity
@property
def fission(self):
return self._fission
@property
def nu_fission(self):
return self._nu_fission
@property
def k_fission(self):
return self._k_fission
@property
def chi(self):
return self._chi
@property
def num_orders(self):
if (self._order is not None) and (self._scatt_type is not None):
if self._scatt_type is 'legendre':
return self._order + 1
else:
return self._order
@name.setter
def name(self, name):
check_type('name for XSdata', name, basestring)
self._name = name
@energy_groups.setter
def energy_groups(self, energy_groups):
# Check validity of energy_groups
check_type("energy_groups", energy_groups, EnergyGroups)
# Check that there is one or more groups
if ((energy_groups.num_groups is None) or
(energy_groups.num_groups < 1)):
msg = 'energy_groups object incorrectly initialized.'
raise ValueError(msg)
self._energy_groups = energy_groups
@representation.setter
def representation(self, representation):
# Check it is of valid type.
check_value('representation', representation, _REPRESENTATIONS)
self._representation = representation
@alias.setter
def alias(self, alias):
if alias is not None:
check_type('alias', alias, basestring)
self._alias = alias
else:
self._alias = self._name
@kT.setter
def kT(self, kT):
# Check validity of type and that the kT value is >= 0
check_type("kT", kT, Real)
check_greater_than("kT", kT, 0.0, equality=True)
self._kT = kT
@scatt_type.setter
def scatt_type(self, scatt_type):
# check to see it is of a valid type and value
check_value("scatt_type", scatt_type, ['legendre', 'histogram',
'tabular'])
self._scatt_type = scatt_type
@order.setter
def order(self, order):
# Check type and value
check_type("order", order, Integral)
check_greater_than("order", order, 0, equality=True)
self._order = order
@tabular_legendre.setter
def tabular_legendre(self, tabular_legendre):
# Check to make sure this is a dict and it has our keys with the
# right values.
check_type("tabular_legendre", tabular_legendre, dict)
if 'enable' in tabular_legendre:
enable = tabular_legendre['enable']
check_type('enable', enable, bool)
else:
msg = "enable must be provided in tabular_legendre"
raise ValueError(msg)
if 'num_points' in tabular_legendre:
num_points = tabular_legendre['num_points']
check_value('num_points', num_points, Integral)
check_greater_than('num_points', num_points, 0)
else:
if not enable:
num_points = 1
else:
num_points = 33
self._tabular_legendre = {'enable': enable, 'num_points': num_points}
@num_polar.setter
def num_polar(self, num_polar):
# Make sure we have positive ints
check_value("num_polar", num_polar, Integral)
check_greater_than("num_polar", num_polar, 0)
self._num_polar = num_polar
@num_azimuthal.setter
def num_azimuthal(self, num_azimuthal):
check_value("num_azimuthal", num_azimuthal, Integral)
check_greater_than("num_azimuthal", num_azimuthal, 0)
self._num_azimuthal = num_azimuthal
@total.setter
def total(self, total):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
# check we have a numpy list
check_type("total", total, np.ndarray, expected_iter_type=Real)
if total.shape == shape:
self._total = np.copy(total)
else:
msg = 'Shape of provided total "{0}" does not match shape ' \
'required, "{1}"'.format(total.shape, shape)
raise ValueError(msg)
@absorption.setter
def absorption(self, absorption):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
# check we have a numpy list
check_type("absorption", absorption, np.ndarray, expected_iter_type=Real)
if absorption.shape == shape:
self._absorption = np.copy(absorption)
else:
msg = 'Shape of provided absorption "{0}" does not match shape ' \
'required, "{1}"'.format(absorption.shape, shape)
raise ValueError(msg)
@fission.setter
def fission(self, fission):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
# check we have a numpy list
check_type("fission", fission, np.ndarray, expected_iter_type=Real)
if fission.shape == shape:
self._fission = np.copy(fission)
if np.sum(self._fission) > 0.0:
self._fissionable = True
else:
msg = 'Shape of provided fission "{0}" does not match shape ' \
'required, "{1}"'.format(fission.shape, shape)
raise ValueError(msg)
@k_fission.setter
def k_fission(self, k_fission):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
# check we have a numpy list
check_type("k_fission", k_fission, np.ndarray, expected_iter_type=Real)
if k_fission.shape == shape:
self._k_fission = np.copy(k_fission)
if np.sum(self._k_fission) > 0.0:
self._fissionable = True
else:
msg = 'Shape of provided k_fission "{0}" does not match shape ' \
'required, "{1}"'.format(k_fission.shape, shape)
raise ValueError(msg)
@chi.setter
def chi(self, chi):
if not self._use_chi:
msg = 'Providing chi when nu_fission already provided as matrix!'
raise ValueError(msg)
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,)
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
# check we have a numpy list
check_type("chi", chi, np.ndarray, expected_iter_type=Real)
if chi.shape == shape:
self._chi = np.copy(chi)
else:
msg = 'Shape of provided chi "{0}" does not match shape ' \
'required, "{1}"'.format(chi.shape, shape)
raise ValueError(msg)
if self._use_chi is not None:
self._use_chi = True
@scatter.setter
def scatter(self, scatter):
if self._representation is 'isotropic':
shape = (self.num_orders, self._energy_groups.num_groups,
self._energy_groups.num_groups)
max_depth = 3
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal, self.num_orders,
self._energy_groups.num_groups,
self._energy_groups.num_groups)
max_depth = 5
# check we have a numpy list
check_iterable_type("scatter", scatter, expected_type=Real,
max_depth=max_depth)
if scatter.shape == shape:
self._scatter = np.copy(scatter)
else:
msg = 'Shape of provided scatter "{0}" does not match shape ' \
'required, "{1}"'.format(scatter.shape, shape)
raise ValueError(msg)
@multiplicity.setter
def multiplicity(self, multiplicity):
if self._representation is 'isotropic':
shape = (self._energy_groups.num_groups,
self._energy_groups.num_groups)
max_depth = 2
elif self._representation is 'angle':
shape = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups,
self._energy_groups.num_groups)
max_depth = 4
# check we have a numpy list
check_iterable_type("multiplicity", multiplicity, expected_type=Real,
max_depth=max_depth)
if multiplicity.shape == shape:
self._multiplicity = np.copy(multiplicity)
else:
msg = 'Shape of provided multiplicity "{0}" does not match shape ' \
'required, "{1}"'.format(multiplicity.shape, shape)
raise ValueError(msg)
@nu_fission.setter
def nu_fission(self, nu_fission):
# nu_fission can be given as a vector or a matrix
# Vector is used when chi also exists.
# Matrix is used when chi does not exist.
# We have to check that the correct form is given, but only if
# chi already has been set. If not, we just check that this is OK
# and set the use_chi flag.
# First lets set our dimensions here since they get used repeatedly
# throughout this code.
if self._representation is 'isotropic':
shape_vec = (self._energy_groups.num_groups,)
shape_mat = (self._energy_groups.num_groups,
self._energy_groups.num_groups)
elif self._representation is 'angle':
shape_vec = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups)
shape_mat = (self._num_polar, self._num_azimuthal,
self._energy_groups.num_groups,
self._energy_groups.num_groups)
# Begin by checking the case when chi has already been given and thus
# the rules for filling in nu_fission are set.
if self._use_chi is not None:
if self._use_chi:
shape = shape_vec
else:
shape = shape_mat
if nu_fission.shape != shape:
msg = "Invalid Shape of Nu_fission!"
raise ValueError(msg)
else:
# Get shape of nu_fission so we can figure if we need chi or not
if nu_fission.shape == shape_vec:
self._use_chi = True
shape = shape_vec
elif nu_fission.shape == shape_mat:
self._use_chi = False
shape = shape_mat
else:
msg = "Invalid Shape of Nu_fission!"
raise ValueError(msg)
# check we have a numpy list
check_type("nu_fission", nu_fission, np.ndarray, expected_iter_type=Real)
self._nu_fission = np.copy(nu_fission)
if np.sum(self._nu_fission) > 0.0:
self._fissionable = True
def _get_xsdata_xml(self):
element = ET.Element("xsdata")
element.set("name", self._name)
if self._alias is not None:
subelement = ET.SubElement(element, 'alias')
subelement.text = self.alias
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)
if self._representation is not None:
subelement = ET.SubElement(element, 'representation')
subelement.text = self._representation
if self._representation == 'angle':
if self._num_azimuthal is not None:
subelement = ET.SubElement(element, 'num_azimuthal')
subelement.text = str(self._num_azimuthal)
if self._num_polar is not None:
subelement = ET.SubElement(element, 'num_polar')
subelement.text = str(self._num_polar)
if self._scatt_type is not None:
subelement = ET.SubElement(element, 'scatt_type')
subelement.text = self._scatt_type
if self._order is not None:
subelement = ET.SubElement(element, 'order')
subelement.text = str(self._order)
if self._tabular_legendre is not None:
subelement = ET.SubElement(element, 'tabular_legendre')
subelement.set('enable', str(self._tabular_legendre['enable']))
subelement.set('num_points', str(self._tabular_legendre['num_points']))
if self._total is not None:
subelement = ET.SubElement(element, 'total')
subelement.text = ndarray_to_string(self._total)
if self._absorption is not None:
subelement = ET.SubElement(element, 'absorption')
subelement.text = ndarray_to_string(self._absorption)
if self._scatter is not None:
subelement = ET.SubElement(element, 'scatter')
subelement.text = ndarray_to_string(self._scatter)
if self._multiplicity is not None:
subelement = ET.SubElement(element, 'multiplicity')
subelement.text = ndarray_to_string(self._multiplicity)
if self._fissionable:
if self._fission is not None:
subelement = ET.SubElement(element, 'fission')
subelement.text = ndarray_to_string(self._fission)
if self._k_fission is not None:
subelement = ET.SubElement(element, 'k_fission')
subelement.text = ndarray_to_string(self._k_fission)
if self._nu_fission is not None:
subelement = ET.SubElement(element, 'nu_fission')
subelement.text = ndarray_to_string(self._nu_fission)
if self._chi is not None:
subelement = ET.SubElement(element, 'chi')
subelement.text = ndarray_to_string(self._chi)
return element
class MGXSLibraryFile(object):
"""Multi-Group Cross Sections file used for an OpenMC simulation.
Corresponds directly to the MG version of the cross_sections.xml input file.
Attributes
----------
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure.
inverse_velocities : Iterable of Real
Inverse of velocities, units of sec/cm
filename : str
XML file to write to.
xsdatas : Iterable of XSdata
Iterable of multi-Group cross section data objects
"""
def __init__(self, energy_groups):
# Initialize MGXSLibraryFile class attributes
self._xsdatas = []
self._energy_groups = energy_groups
self._inverse_velocities = None
self._cross_sections_file = ET.Element("cross_sections")
@property
def inverse_velocities(self):
return self._inverse_velocities
@property
def energy_groups(self):
return self._energy_groups
@inverse_velocities.setter
def inverse_velocities(self, inverse_velocities):
cv.check_type('inverse_velocities', inverse_velocities, Iterable, Real)
cv.check_greater_than('number of inverse_velocities',
len(inverse_velocities), 0.0)
self._inverse_velocities = np.array(inverse_velocities)
@energy_groups.setter
def energy_groups(self, energy_groups):
check_type("energy groups", energy_groups, EnergyGroups)
self._energy_groups = energy_groups
def add_xsdata(self, xsdata):
"""Add an XSdata entry to the file.
Parameters
----------
xsdata : XSdata
MGXS information to add
"""
# Check the type
if not isinstance(xsdata, XSdata):
msg = 'Unable to add a non-XSdata "{0}" to the ' \
'MGXSLibraryFile'.format(xsdata)
raise ValueError(msg)
# Make sure energy groups match.
if xsdata.energy_groups != self._energy_groups:
msg = 'Energy groups of XSdata do not match that of MGXSLibraryFile!'
raise ValueError(msg)
self._xsdatas.append(xsdata)
def add_xsdatas(self, xsdatas):
"""Add multiple xsdatas to the file.
Parameters
----------
xsdatas : tuple or list of XSdata
XSdatas to add
"""
if not isinstance(xsdatas, Iterable):
msg = 'Unable to create OpenMC xsdatas.xml file from "{0}" which ' \
'is not iterable'.format(xsdatas)
raise ValueError(msg)
for xsdata in xsdatas:
self.add_xsdata(xsdata)
def remove_xsdata(self, xsdata):
"""Remove a xsdata from the file
Parameters
----------
xsdata : XSdata
XSdata to remove
"""
if not isinstance(xsdata, XSdata):
msg = 'Unable to remove a non-XSdata "{0}" from the ' \
'XSdatasFile'.format(xsdata)
raise ValueError(msg)
self._xsdatas.remove(xsdata)
def _create_groups_subelement(self):
if self._energy_groups is not None:
element = ET.SubElement(self._cross_sections_file, "groups")
element.text = str(self._energy_groups.num_groups)
def _create_group_structure_subelement(self):
if self._energy_groups is not None:
element = ET.SubElement(self._cross_sections_file,
"group_structure")
element.text = ' '.join(map(str, self._energy_groups.group_edges))
def _create_inverse_velocities_subelement(self):
if self._inverse_velocities is not None:
element = ET.SubElement(self._cross_sections_file,
"inverse_velocities")
element.text = ' '.join(map(str, self._inverse_velocities))
def _create_xsdata_subelements(self):
for xsdata in self._xsdatas:
xml_element = xsdata._get_xsdata_xml()
self._cross_sections_file.append(xml_element)
def export_to_xml(self, filename='mg_cross_sections.xml'):
"""Create an mg_cross_sections.xml file that can be used for a
simulation.
Parameters
----------
filename : str, optional
filename of file, default is mg_cross_sections.xml
"""
# Reset xml element tree
self._cross_sections_file.clear()
self._create_groups_subelement()
self._create_group_structure_subelement()
self._create_inverse_velocities_subelement()
self._create_xsdata_subelements()
# Clean the indentation in the file to be user-readable
sort_xml_elements(self._cross_sections_file)
clean_xml_indentation(self._cross_sections_file)
# Write the XML Tree to the xsdatas.xml file
tree = ET.ElementTree(self._cross_sections_file)
tree.write(filename, xml_declaration=True,
encoding='utf-8', method="xml")

50
openmc/settings.py
View file

@ -70,11 +70,15 @@ class SettingsFile(object):
cross_sections : str
Indicates the path to an XML cross section listing file (usually named
cross_sections.xml). If it is not set, the :envvar:`CROSS_SECTIONS`
environment variable will be used to find the path to the XML cross
section listing.
energy_grid : str
Set the method used to search energy grids. Acceptable values are
'nuclide', 'logarithm', and 'material-union'.
environment variable will be used for continuous-energy calculations
and :envvar:`MG_CROSS_SECTIONS` will be used for multi-group
calculations to find the path to the XML cross section file.
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
Maximum scattering order to apply globally when in multi-group mode.
ptables : bool
Determine whether probability tables are used.
run_cmfd : bool
@ -134,6 +138,10 @@ class SettingsFile(object):
self._particles = None
self._keff_trigger = None
# Energy mode subelement
self._energy_mode = None
self._max_order = None
# Source subelement
self._source = None
@ -221,6 +229,14 @@ class SettingsFile(object):
def keff_trigger(self):
return self._keff_trigger
@property
def energy_mode(self):
return self._energy_mode
@property
def max_order(self):
return self._max_order
@property
def source(self):
return self._source
@ -438,6 +454,18 @@ class SettingsFile(object):
self._keff_trigger = keff_trigger
@energy_mode.setter
def energy_mode(self, energy_mode):
check_value('energy mode', energy_mode,
['continuous-energy', 'multi-group'])
self._energy_mode = energy_mode
@max_order.setter
def max_order(self, max_order):
check_type('maximum scattering order', max_order, Integral)
check_greater_than('maximum scattering order', max_order, 0, True)
self._max_order = max_order
@source.setter
def source(self, source):
if isinstance(source, Source):
@ -782,6 +810,16 @@ class SettingsFile(object):
subelement = ET.SubElement(element, key)
subelement.text = str(self._keff_trigger[key]).lower()
def _create_energy_mode_subelement(self):
if self._energy_mode is not None:
element = ET.SubElement(self._settings_file, "energy_mode")
element.text = str(self._energy_mode)
def _create_max_order_subelement(self):
if self._max_order is not None:
element = ET.SubElement(self._settings_file, "max_order")
element.text = str(self._max_order)
def _create_source_subelement(self):
if self.source is not None:
for source in self.source:
@ -1025,6 +1063,8 @@ class SettingsFile(object):
self._create_confidence_intervals()
self._create_cross_sections_subelement()
self._create_energy_grid_subelement()
self._create_energy_mode_subelement()
self._create_max_order_subelement()
self._create_ptables_subelement()
self._create_run_cmfd_subelement()
self._create_seed_subelement()

4
openmc/statepoint.py
View file

@ -103,11 +103,11 @@ class StatePoint(object):
raise IOError('Could not read statepoint file. This most likely '
'means the statepoint file was produced by a different '
'version of OpenMC than the one you are using.')
if self._f['revision'].value != 14:
if self._f['revision'].value != 15:
raise IOError('Statepoint file has a file revision of {} '
'which is not consistent with the revision this '
'version of OpenMC expects ({}).'.format(
self._f['revision'].value, 14))
self._f['revision'].value, 15))
# Set flags for what data has been read
self._meshes_read = False

5
openmc/summary.py
View file

@ -60,6 +60,9 @@ class Summary(object):
# Read date and time
self.date_and_time = self._f['date_and_time'][...]
# Read if continuous-energy or multi-group
self.run_CE = (self._f['run_CE'].value == 1)
self.n_batches = self._f['n_batches'].value
self.n_particles = self._f['n_particles'].value
self.n_active = self._f['n_active'].value
@ -276,7 +279,7 @@ class Summary(object):
# Get the distribcell index
ind = self._f['geometry/cells'][key]['distribcell_index'].value
if ind != 0:
cell.distribcell_index = ind
cell.distribcell_index = ind
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell

32
src/ace.F90
View file

@ -1,6 +1,6 @@
module ace
use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing
use ace_header, only: Reaction
use constants
use distribution_univariate, only: Uniform, Equiprobable, Tabular
use endf, only: is_fission, is_disappearance
@ -11,7 +11,9 @@ module ace
use global
use list_header, only: ListInt
use material_header, only: Material
use nuclide_header
use output, only: write_message
use sab_header
use set_header, only: SetChar
use secondary_header, only: AngleEnergy
use secondary_correlated, only: CorrelatedAngleEnergy
@ -33,11 +35,11 @@ module ace
contains
!===============================================================================
! READ_XS reads all the cross sections for the problem and stores them in
! READ_ACE_XS reads all the cross sections for the problem and stores them in
! nuclides and sab_tables arrays
!===============================================================================
subroutine read_xs()
subroutine read_ace_xs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
@ -52,7 +54,7 @@ contains
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
type(Material), pointer :: mat
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(SetChar) :: already_read
@ -231,7 +233,7 @@ contains
end if
end do
end subroutine read_xs
end subroutine read_ace_xs
!===============================================================================
! READ_ACE_TABLE reads a single cross section table in either ASCII or binary
@ -263,9 +265,9 @@ contains
character(10) :: mat ! material identifier
character(70) :: comment ! comment for ACE table
character(MAX_FILE_LEN) :: filename ! path to ACE cross section library
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(XsListing), pointer :: listing
type(XsListing), pointer :: listing
! determine path, record length, and location of table
listing => xs_listings(i_listing)
@ -420,8 +422,8 @@ contains
!===============================================================================
subroutine read_esz(nuc, data_0K)
type(Nuclide), intent(inout) :: nuc
logical, intent(in) :: data_0K ! are we reading 0K data?
type(NuclideCE), intent(inout) :: nuc
logical, intent(in) :: data_0K ! are we reading 0K data?
integer :: NE ! number of energy points for total and elastic cross sections
integer :: i ! index in 0K elastic xs array for this nuclide
@ -508,7 +510,7 @@ contains
!===============================================================================
subroutine read_nu_data(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: i ! loop index
integer :: JXS2 ! location for fission nu data
@ -712,7 +714,7 @@ contains
!===============================================================================
subroutine read_reactions(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: i ! loop indices
integer :: i_fission ! index in nuc % index_fission
@ -892,7 +894,7 @@ contains
!===============================================================================
subroutine read_angular_dist(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: LOCB ! location of angular distribution for given MT
integer :: NE ! number of incoming energies
@ -996,7 +998,7 @@ contains
!===============================================================================
subroutine read_energy_dist(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: i ! loop index
integer :: n
@ -1384,7 +1386,7 @@ contains
!===============================================================================
subroutine read_unr_res(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: JXS23 ! location of URR data
integer :: lc ! locator
@ -1472,7 +1474,7 @@ contains
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide), intent(inout) :: nuc
type(NuclideCE), intent(inout) :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy

230
src/ace_header.F90
View file

@ -24,9 +24,8 @@ module ace_header
real(8), allocatable :: sigma(:) ! Cross section values
type(SecondaryDistribution) :: secondary
! Type-Bound procedures
contains
procedure :: clear => reaction_clear ! Deallocates Reaction
contains
procedure :: clear => reaction_clear ! Deallocates Reaction
end type Reaction
!===============================================================================
@ -44,210 +43,6 @@ module ace_header
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
end type UrrData
!===============================================================================
! NUCLIDE contains all the data for an ACE-format continuous-energy cross
! section. The ACE format (A Compact ENDF format) is used in MCNP and several
! other Monte Carlo codes.
!===============================================================================
type Nuclide
character(10) :: name ! name of nuclide, e.g. 92235.03c
integer :: zaid ! Z and A identifier, e.g. 92235
integer :: listing ! index in xs_listings
real(8) :: awr ! weight of nucleus in neutron masses
real(8) :: kT ! temperature in MeV (k*T)
! Linked list of indices in nuclides array of instances of this same nuclide
type(VectorInt) :: nuc_list
! Energy grid information
integer :: n_grid ! # of nuclide grid points
integer, allocatable :: grid_index(:) ! log grid mapping indices
real(8), allocatable :: energy(:) ! energy values corresponding to xs
! Microscopic cross sections
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
! Resonance scattering info
logical :: resonant = .false. ! resonant scatterer?
character(10) :: name_0K = '' ! name of 0K nuclide, e.g. 92235.00c
character(16) :: scheme ! target velocity sampling scheme
integer :: n_grid_0K ! number of 0K energy grid points
real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs
real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
real(8) :: E_min ! lower cutoff energy for res scattering
real(8) :: E_max ! upper cutoff energy for res scattering
! Fission information
logical :: fissionable ! nuclide is fissionable?
logical :: has_partial_fission ! nuclide has partial fission reactions?
integer :: n_fission ! # of fission reactions
integer, allocatable :: index_fission(:) ! indices in reactions
! Total fission neutron emission
integer :: nu_t_type
real(8), allocatable :: nu_t_data(:)
! Prompt fission neutron emission
integer :: nu_p_type
real(8), allocatable :: nu_p_data(:)
! Delayed fission neutron emission
integer :: nu_d_type
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
! Unresolved resonance data
logical :: urr_present
integer :: urr_inelastic
type(UrrData), pointer :: urr_data => null()
! Reactions
integer :: n_reaction ! # of reactions
type(Reaction), allocatable :: reactions(:)
type(DictIntInt) :: reaction_index ! map MT values to index in reactions
! array; used at tally-time
! Type-Bound procedures
contains
procedure :: clear => nuclide_clear ! Deallocates Nuclide
end type Nuclide
!===============================================================================
! NUCLIDE0K temporarily contains all 0K cross section data and other parameters
! needed to treat resonance scattering before transferring them to NUCLIDE
!===============================================================================
type Nuclide0K
character(10) :: nuclide ! name of nuclide, e.g. U-238
character(16) :: scheme = 'ares' ! target velocity sampling scheme
character(10) :: name ! name of nuclide, e.g. 92235.03c
character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c
real(8) :: E_min = 0.01e-6_8 ! lower cutoff energy for res scattering
real(8) :: E_max = 1000.0e-6_8 ! upper cutoff energy for res scattering
end type Nuclide0K
!===============================================================================
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
! thermal scattering collisions which utilize a continuous secondary energy
! representation.
!===============================================================================
type DistEnergySab
integer :: n_e_out
real(8), allocatable :: e_out(:)
real(8), allocatable :: e_out_pdf(:)
real(8), allocatable :: e_out_cdf(:)
real(8), allocatable :: mu(:,:)
end type DistEnergySab
!===============================================================================
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
! of light isotopes such as water, graphite, Be, etc
!===============================================================================
type SAlphaBeta
character(10) :: 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
! threshold for S(a,b) treatment (usually ~4 eV)
real(8) :: threshold_inelastic
real(8) :: threshold_elastic = ZERO
! Inelastic scattering data
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
real(8), allocatable :: inelastic_e_out(:,:)
real(8), allocatable :: inelastic_mu(:,:,:)
! The following is used only if secondary_mode is 3
! The different implementation is necessary because the continuous
! representation has a variable number of outgoing energy points for each
! incoming energy
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
! Elastic scattering data
integer :: elastic_mode ! elastic mode (discrete/exact)
integer :: n_elastic_e_in ! # of incoming E for elastic
integer :: n_elastic_mu ! # of outgoing angles for elastic
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
real(8), allocatable :: elastic_mu(:,:)
end type SAlphaBeta
!===============================================================================
! XSLISTING contains data read from a cross_sections.xml file
!===============================================================================
type XsListing
character(12) :: name ! table name, e.g. 92235.70c
character(12) :: alias ! table alias, e.g. U-235.70c
integer :: type ! type of table (cont-E neutron, S(A,b), etc)
integer :: zaid ! ZAID identifier = 1000*Z + A
integer :: filetype ! ASCII or BINARY
integer :: location ! location of table within library
integer :: recl ! record length for library
integer :: entries ! number of entries per record
real(8) :: awr ! atomic weight ratio (# of neutron masses)
real(8) :: kT ! Boltzmann constant * temperature (MeV)
logical :: metastable ! is this nuclide metastable?
character(MAX_FILE_LEN) :: path ! path to library containing table
end type XsListing
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
!===============================================================================
type NuclideMicroXS
integer :: index_grid ! index on nuclide energy grid
integer :: index_temp ! temperature index for nuclide
real(8) :: last_E = ZERO ! last evaluated energy
real(8) :: interp_factor ! interpolation factor on nuc. energy grid
real(8) :: total ! microscropic total xs
real(8) :: elastic ! microscopic elastic scattering xs
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
! 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
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
real(8) :: last_prn
end type NuclideMicroXS
!===============================================================================
! MATERIALMACROXS contains cached macroscopic cross sections for the material a
! particle is traveling through
!===============================================================================
type MaterialMacroXS
real(8) :: total ! macroscopic total xs
real(8) :: elastic ! macroscopic elastic scattering xs
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
end type MaterialMacroXS
contains
!===============================================================================
@ -260,25 +55,4 @@ module ace_header
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
end subroutine reaction_clear
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide.
!===============================================================================
subroutine nuclide_clear(this)
class(Nuclide), intent(inout) :: this
integer :: i ! Loop counter
if (associated(this % urr_data)) deallocate(this % urr_data)
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()
end do
end if
call this % reaction_index % clear()
end subroutine nuclide_clear
end module ace_header

2
src/bank_header.F90
View file

@ -14,7 +14,7 @@ module bank_header
real(C_DOUBLE) :: wgt ! weight of bank site
real(C_DOUBLE) :: xyz(3) ! location of bank particle
real(C_DOUBLE) :: uvw(3) ! diretional cosines
real(C_DOUBLE) :: E ! energy
real(C_DOUBLE) :: E ! energy / energy group if in MG mode.
integer(C_INT) :: delayed_group ! delayed group
end type Bank

18
src/cmfd_data.F90
View file

@ -49,13 +49,13 @@ contains
subroutine compute_xs()
use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, &
FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, &
ONE, TINY_BIT
use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, &
FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, &
OUT_FRONT, IN_TOP, OUT_TOP, CMFD_NOACCEL, ZERO, &
ONE, TINY_BIT
use error, only: fatal_error
use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes,&
matching_bins
matching_bins
use mesh, only: mesh_indices_to_bin
use mesh_header, only: RegularMesh
use string, only: to_str
@ -625,10 +625,10 @@ contains
subroutine compute_dhat()
use constants, only: CMFD_NOACCEL, ZERO
use global, only: cmfd, cmfd_coremap, dhat_reset
use output, only: write_message
use string, only: to_str
use constants, only: CMFD_NOACCEL, ZERO
use global, only: cmfd, cmfd_coremap, dhat_reset
use output, only: write_message
use string, only: to_str
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction

6
src/cmfd_execute.F90
View file

@ -89,9 +89,9 @@ contains
subroutine calc_fission_source()
use constants, only: CMFD_NOACCEL, ZERO, TWO
use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch
use string, only: to_str
use constants, only: CMFD_NOACCEL, ZERO, TWO
use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch
use string, only: to_str
#ifdef MPI
use global, only: mpi_err

25
src/cmfd_input.F90
View file

@ -45,14 +45,14 @@ contains
subroutine read_cmfd_xml()
use constants, only: ZERO, ONE
use error, only: fatal_error, warning
use error, only: fatal_error, warning
use global
use output, only: write_message
use string, only: to_lower
use output, only: write_message
use string, only: to_lower
use xml_interface
use, intrinsic :: ISO_FORTRAN_ENV
integer :: i
integer :: i, g
integer :: ng
integer :: n_params
integer, allocatable :: iarray(:)
@ -102,6 +102,23 @@ contains
if(.not.allocated(cmfd%egrid)) allocate(cmfd%egrid(ng))
call get_node_array(node_mesh, "energy", cmfd%egrid)
cmfd % indices(4) = ng - 1 ! sets energy group dimension
! If using MG mode, check to see if these egrid points at least match
! the MG Data breakpoints
if (.not. run_CE) then
do i = 1, ng
found = .false.
do g = 1, energy_groups + 1
if (cmfd % egrid(i) == energy_bins(g)) then
found = .true.
exit
end if
end do
if (.not. found) then
call fatal_error("CMFD energy mesh boundaries must align with&
& boundaries of multi-group data!")
end if
end do
end if
else
if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
cmfd % egrid = [ ZERO, 20.0_8 ]

20
src/constants.F90
View file

@ -11,10 +11,10 @@ module constants
integer, parameter :: VERSION_RELEASE = 1
! Revision numbers for binary files
integer, parameter :: REVISION_STATEPOINT = 14
integer, parameter :: REVISION_STATEPOINT = 15
integer, parameter :: REVISION_PARTICLE_RESTART = 1
integer, parameter :: REVISION_TRACK = 1
integer, parameter :: REVISION_SUMMARY = 2
integer, parameter :: REVISION_SUMMARY = 3
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -161,9 +161,14 @@ module constants
! Angular distribution type
integer, parameter :: &
ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
ANGLE_32_EQUI = 2, & ! 32 equiprobable bins
ANGLE_TABULAR = 3 ! Tabular angular distribution
ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution (CE)
ANGLE_32_EQUI = 2, & ! 32 equiprobable bins (CE)
ANGLE_TABULAR = 3, & ! Tabular angular distribution (CE or MG)
ANGLE_LEGENDRE = 4, & ! Legendre angular distribution (MG)
ANGLE_HISTOGRAM = 5 ! Histogram angular distribution (MG)
! Number of mu bins to use when converting Legendres to tabular type
integer, parameter :: DEFAULT_NMU = 33
! Secondary energy mode for S(a,b) inelastic scattering
integer, parameter :: &
@ -207,6 +212,11 @@ module constants
ACE_THERMAL = 2, & ! thermal S(a,b) scattering data
ACE_DOSIMETRY = 3 ! dosimetry cross sections
! MGXS Table Types
integer, parameter :: &
MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data
MGXS_ANGLE = 2 ! Data by Angular Bins
! Fission neutron emission (nu) type
integer, parameter :: &
NU_NONE = 0, & ! No nu values (non-fissionable)

18
src/cross_section.F90
View file

@ -1,6 +1,6 @@
module cross_section
use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData
use ace_header, only: Reaction, UrrData
use constants
use energy_grid, only: grid_method, log_spacing
use error, only: fatal_error
@ -8,8 +8,10 @@ module cross_section
use global
use list_header, only: ListElemInt
use material_header, only: Material
use nuclide_header
use particle_header, only: Particle
use random_lcg, only: prn
use sab_header, only: SAlphaBeta
use search, only: binary_search
implicit none
@ -146,8 +148,8 @@ contains
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
type(NuclideCE), pointer :: nuc
type(Material), pointer :: mat
! Set pointer to nuclide and material
nuc => nuclides(i_nuclide)
@ -364,8 +366,8 @@ contains
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer :: urr
type(Nuclide), pointer :: nuc
type(UrrData), pointer :: urr
type(NuclideCE), pointer :: nuc
micro_xs(i_nuclide) % use_ptable = .true.
@ -531,9 +533,9 @@ contains
!===============================================================================
pure function elastic_xs_0K(E, nuc) result(xs_out)
real(8), intent(in) :: E ! trial energy
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
real(8), intent(in) :: E ! trial energy
type(NuclideCE), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid

6
src/distribution_multivariate.F90
View file

@ -1,9 +1,9 @@
module distribution_multivariate
use constants, only: ONE, TWO, PI
use constants, only: ONE, TWO, PI
use distribution_univariate, only: Distribution
use math, only: rotate_angle
use random_lcg, only: prn
use random_lcg, only: prn
use math, only: rotate_angle
implicit none

8
src/distribution_univariate.F90
View file

@ -1,11 +1,11 @@
module distribution_univariate
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use error, only: fatal_error
use random_lcg, only: prn
use string, only: to_lower
use math, only: maxwell_spectrum, watt_spectrum
use string, only: to_lower
use xml_interface
implicit none

17
src/eigenvalue.F90
View file

@ -4,16 +4,15 @@ module eigenvalue
use message_passing
#endif
use constants, only: ZERO
use error, only: fatal_error, warning
use constants, only: ZERO
use error, only: fatal_error, warning
use global
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: RegularMesh
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
implicit none

10
src/energy_distribution.F90
View file

@ -1,11 +1,11 @@
module energy_distribution
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tab1
use interpolation, only: interpolate_tab1
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
use search, only: binary_search
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

6
src/energy_grid.F90
View file

@ -27,7 +27,7 @@ contains
integer :: i ! index in nuclides array
integer :: j ! index in materials array
type(ListReal) :: list
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
type(Material), pointer :: mat
call write_message("Creating unionized energy grid...", 5)
@ -70,7 +70,7 @@ contains
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
type(NuclideCE), pointer :: nuc
! Set minimum/maximum energies
E_max = energy_max_neutron
@ -179,7 +179,7 @@ contains
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(NuclideCE), pointer :: nuc
type(Material), pointer :: mat
do k = 1, n_materials

16
src/fission.F90
View file

@ -1,10 +1,10 @@
module fission
use ace_header, only: Nuclide
use nuclide_header, only: NuclideCE
use constants
use error, only: fatal_error
use interpolation, only: interpolate_tab1
use search, only: binary_search
use error, only: fatal_error
use interpolation, only: interpolate_tab1
use search, only: binary_search
implicit none
@ -16,7 +16,7 @@ contains
!===============================================================================
pure function nu_total(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
@ -49,7 +49,7 @@ contains
!===============================================================================
pure function nu_prompt(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
@ -86,7 +86,7 @@ contains
!===============================================================================
pure function nu_delayed(nuc, E) result(nu)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
@ -109,7 +109,7 @@ contains
!===============================================================================
pure function yield_delayed(nuc, E, g) result(yield)
type(Nuclide), intent(in) :: nuc ! nuclide from which to find nu
type(NuclideCE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: yield ! delayed neutron precursor yield
integer, intent(in) :: g ! the delayed neutron precursor group

69
src/global.F90
View file

@ -1,15 +1,16 @@
module global
use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, &
MaterialMacroXS, Nuclide0K
use bank_header, only: Bank
use cmfd_header
use constants
use dict_header, only: DictCharInt, DictIntInt
use geometry_header, only: Cell, Universe, Lattice, LatticeContainer
use macroxs_header, only: MacroXSContainer
use material_header, only: Material
use mesh_header, only: RegularMesh
use nuclide_header
use plot_header, only: ObjectPlot
use sab_header, only: SAlphaBeta
use set_header, only: SetInt
use surface_header, only: SurfaceContainer
use source_header, only: SourceDistribution
@ -58,39 +59,75 @@ module global
integer :: n_lost_particles
! ============================================================================
! CROSS SECTION RELATED VARIABLES
! ENERGY TREATMENT RELATED VARIABLES
logical :: run_CE = .true. ! Run in CE mode?
! ============================================================================
! CROSS SECTION RELATED VARIABLES NEEDED REGARDLESS OF CE OR MG
! Cross section arrays
type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections
type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables
type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings
integer :: n_nuclides_total ! Number of nuclide cross section tables
integer :: n_listings ! Number of listings in cross_sections.xml
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
type(MaterialMacroXS) :: material_xs ! Cache for current material
integer :: n_nuclides_total ! Number of nuclide cross section tables
! Dictionaries to look up cross sections and listings
type(DictCharInt) :: nuclide_dict
type(DictCharInt) :: xs_listing_dict
! Default xs identifier (e.g. 70c or 300K)
character(5):: default_xs
! ============================================================================
! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES
! Cross section arrays
type(NuclideCE), allocatable, target :: nuclides(:) ! Nuclide cross-sections
type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables
integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables
integer :: n_listings ! Number of listings in cross_sections.xml
! Minimum/maximum energies
real(8) :: energy_min_neutron = ZERO
real(8) :: energy_max_neutron = INFINITY
! Dictionaries to look up cross sections and listings
type(DictCharInt) :: nuclide_dict
type(DictCharInt) :: sab_dict
type(DictCharInt) :: xs_listing_dict
! Unreoslved resonance probablity tables
logical :: urr_ptables_on = .true.
! Default xs identifier (e.g. 70c)
character(3):: default_xs
! What to assume for expanding natural elements
integer :: default_expand = ENDF_BVII1
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
! Cross section arrays
type(NuclideMGContainer), allocatable, target :: nuclides_MG(:)
! Cross section caches
type(MacroXSContainer), target, allocatable :: macro_xs(:)
! Number of energy groups
integer :: energy_groups
! Energy group structure
real(8), allocatable :: energy_bins(:)
! Midpoint of the energy group structure
real(8), allocatable :: energy_bin_avg(:)
! Inverse velocities of the energy groups (provided or estimated)
real(8), allocatable :: inverse_velocities(:)
! Maximum Data Order
integer :: max_order
! ============================================================================
! TALLY-RELATED VARIABLES
@ -444,6 +481,14 @@ contains
deallocate(nuclides_0K)
end if
if (allocated(nuclides_MG)) then
deallocate(nuclides_MG)
end if
if (allocated(macro_xs)) then
deallocate(macro_xs)
end if
if (allocated(sab_tables)) deallocate(sab_tables)
if (allocated(xs_listings)) deallocate(xs_listings)
if (allocated(micro_xs)) deallocate(micro_xs)

90
src/initialize.F90
View file

@ -1,30 +1,31 @@
module initialize
use ace, only: read_xs, same_nuclide_list
use bank_header, only: Bank
use ace, only: read_ace_xs, same_nuclide_list
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method, unionized_grid
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
&BASE_UNIVERSE
use global
use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,&
hdf5_tallyresult_t, hdf5_integer8_t
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use output, only: title, header, print_version, write_message, &
print_usage, write_xs_summary, print_plot
use random_lcg, only: initialize_prng
use state_point, only: load_state_point
use string, only: to_str, str_to_int, starts_with, ends_with
use summary, only: write_summary
use tally_header, only: TallyObject, TallyResult, TallyFilter
use tally_initialize, only: configure_tallies
use hdf5_interface, only: file_open, read_dataset, file_close, hdf5_bank_t,&
hdf5_tallyresult_t, hdf5_integer8_t
use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml
use material_header, only: Material
use mgxs_data, only: read_mgxs, same_NuclideMG_list, create_macro_xs
use output, only: title, header, print_version, write_message, &
print_usage, write_xs_summary, 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
use summary, only: write_summary
use tally_header, only: TallyObject, TallyResult, TallyFilter
use tally_initialize,only: configure_tallies
use tally, only: init_tally_routines
#ifdef MPI
use message_passing
@ -114,27 +115,44 @@ contains
! Read ACE-format cross sections
call time_read_xs%start()
call read_xs()
if (run_CE) then
call read_ace_xs()
else
call read_mgxs()
end if
call time_read_xs%stop()
! Create linked lists for multiple instances of the same nuclide
call same_nuclide_list()
if (run_CE) then
call same_nuclide_list()
else
call same_nuclidemg_list()
end if
! 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 information needed for nuclear data
if (run_CE) then
! 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
else
! Create material macroscopic data for MGXS
call create_macro_xs()
end if
! Allocate and setup tally stride, matching_bins, and tally maps
call configure_tallies()
! Set up tally procedure pointers
call init_tally_routines()
! Determine how much work each processor should do
call calculate_work()

635
src/input_xml.F90
View file

File diff suppressed because it is too large Load diff

860
src/macroxs_header.F90 Normal file
View file

@ -0,0 +1,860 @@
module macroxs_header
use constants, only: MAX_FILE_LEN, ZERO, ONE, TWO, PI
use list_header, only: ListInt
use material_header, only: material
use math, only: calc_pn, calc_rn, expand_harmonic, find_angle
use nuclide_header
use random_lcg, only: prn
use scattdata_header
implicit none
!===============================================================================
! MACROXS_* contains cached macroscopic cross sections for the material a
! particle is traveling through
!===============================================================================
type, abstract :: MacroXS
! Data Order
integer :: order
contains
procedure(macroxs_init_), deferred :: init ! initializes object
procedure(macroxs_get_xs_), deferred :: get_xs ! Return xs
! Sample the outgoing energy from a fission event
procedure(macroxs_sample_fission_), deferred :: sample_fission_energy
! Sample the outgoing energy and angle from a scatter event
procedure(macroxs_sample_scatter_), deferred :: sample_scatter
! Calculate the material specific MGXS data from the nuclides
procedure(macroxs_calculate_xs_), deferred :: calculate_xs
end type MacroXS
abstract interface
subroutine macroxs_init_(this, mat, nuclides, groups, get_kfiss, get_fiss, &
max_order, scatt_type, legendre_mu_points, &
error_code, error_text)
import MacroXS, Material, NuclideMGContainer, MAX_LINE_LEN
class(MacroXS), intent(inout) :: this ! The MacroXS to initialize
type(Material), pointer, intent(in) :: mat ! base material
type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
logical, intent(in) :: get_kfiss ! Should we get kfiss data?
logical, intent(in) :: get_fiss ! Should we get fiss data?
integer, intent(in) :: max_order ! Maximum requested order
integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt?
integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular?
integer, intent(inout) :: error_code ! Code signifying error
character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print
end subroutine macroxs_init_
function macroxs_get_xs_(this, g, xstype, gout, uvw) result(xs)
import MacroXS
class(MacroXS), intent(in) :: this ! The MacroXS to initialize
integer, intent(in) :: g ! Incoming Energy group
character(*) , intent(in) :: xstype ! Cross Section Type
integer, optional, intent(in) :: gout ! Outgoing Energy group
real(8), optional, intent(in) :: uvw(3) ! Requested Angle
real(8) :: xs ! Resultant xs
end function macroxs_get_xs_
function macroxs_sample_fission_(this, gin, uvw) result(gout)
import MacroXS
class(MacroXS), intent(in) :: this ! Data to work with
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
end function macroxs_sample_fission_
subroutine macroxs_sample_scatter_(this, uvw, gin, gout, mu, wgt)
import MacroXS
class(MacroXS), intent(in) :: this
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
end subroutine macroxs_sample_scatter_
subroutine macroxs_calculate_xs_(this, gin, uvw, xs)
import MacroXS, MaterialMacroXS
class(MacroXS), intent(in) :: this
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs
end subroutine macroxs_calculate_xs_
end interface
type, extends(MacroXS) :: MacroXSIso
! Microscopic cross sections
real(8), allocatable :: total(:) ! total cross section
real(8), allocatable :: absorption(:) ! absorption cross section
class(ScattData), allocatable :: scatter ! scattering information
real(8), allocatable :: nu_fission(:) ! nu-fission
real(8), allocatable :: k_fission(:) ! kappa-fission
real(8), allocatable :: fission(:) ! fission x/s
real(8), allocatable :: scattxs(:) ! scattering xs
real(8), allocatable :: chi(:,:) ! fission spectra
contains
procedure :: init => macroxsiso_init ! inits object
procedure :: get_xs => macroxsiso_get_xs ! Returns xs
procedure :: sample_fission_energy => macroxsiso_sample_fission_energy
procedure :: sample_scatter => macroxsiso_sample_scatter
procedure :: calculate_xs => macroxsiso_calculate_xs
end type MacroXSIso
type, extends(MacroXS) :: MacroXSAngle
! Macroscopic cross sections
real(8), allocatable :: total(:,:,:) ! total cross section
real(8), allocatable :: absorption(:,:,:) ! absorption cross section
type(ScattDataContainer), allocatable :: scatter(:,:) ! scattering information
real(8), allocatable :: nu_fission(:,:,:) ! nu-fission
real(8), allocatable :: k_fission(:,:,:) ! kappa-fission
real(8), allocatable :: fission(:,:,:) ! fission x/s
real(8), allocatable :: chi(:,:,:,:) ! fission spectra
real(8), allocatable :: scattxs(:,:,:) ! scattering xs
real(8), allocatable :: polar(:) ! polar angles
real(8), allocatable :: azimuthal(:) ! azimuthal angles
contains
procedure :: init => macroxsangle_init ! inits object
procedure :: get_xs => macroxsangle_get_xs ! Returns xs
procedure :: sample_fission_energy => macroxsangle_sample_fission_energy
procedure :: sample_scatter => macroxsangle_sample_scatter
procedure :: calculate_xs => macroxsangle_calculate_xs
end type MacroXSAngle
!===============================================================================
! MACROXSCONTAINER pointer array for storing MacroXS objects.
!===============================================================================
type MacroXSContainer
class(MacroXS), allocatable :: obj
end type MacroXSContainer
contains
!===============================================================================
! MACROXS*_INIT sets the MacroXS Data
!===============================================================================
subroutine macroxsiso_init(this, mat, nuclides, groups, get_kfiss, get_fiss, &
max_order, scatt_type, legendre_mu_points, error_code, error_text)
class(MacroXSIso), intent(inout) :: this ! The MacroXS to initialize
type(Material), pointer, intent(in) :: mat ! base material
type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
logical, intent(in) :: get_kfiss ! Should we get kfiss data?
logical, intent(in) :: get_fiss ! Should we get fiss data?
integer, intent(in) :: max_order ! Maximum requested order
integer, intent(in) :: scatt_type ! How is data presented
integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular?
integer, intent(inout) :: error_code ! Code signifying error
character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print
integer :: i ! loop index over nuclides
integer :: gin, gout ! group indices
real(8) :: atom_density ! atom density of a nuclide
integer :: imu
real(8) :: norm
integer :: mat_max_order, order, l
real(8), allocatable :: temp_mult(:,:)
real(8), allocatable :: temp_energy(:,:)
real(8), allocatable :: scatt_coeffs(:,:,:)
! Initialize error data
error_code = 0
error_text = ''
! If we have tabular only data, then make sure all datasets have same size
if (scatt_type == ANGLE_HISTOGRAM) then
! Check all scattering data of same size
order = nuclides(mat % nuclide(1)) % obj % order
do i = 2, mat % n_nuclides
if (order /= nuclides(mat % nuclide(i)) % obj % order) then
error_code = 1
error_text = "All Histogram Scattering Entries Must Be Same Length!"
return
end if
end do
! Ok, got our order, store it
this % order = order
! Allocate stuff for later
allocate(scatt_coeffs(order, groups, groups))
scatt_coeffs = ZERO
allocate(ScattDataHistogram :: this % scatter)
else if (scatt_type == ANGLE_TABULAR) then
! Check all scattering data of same size
order = nuclides(mat % nuclide(1)) % obj % order
do i = 2, mat % n_nuclides
if (order /= nuclides(mat % nuclide(i)) % obj % order) then
error_code = 1
error_text = "All Tabular Scattering Entries Must Be Same Length!"
return
end if
end do
! Ok, got our order, store it
this % order = order
! Allocate stuff for later
allocate(scatt_coeffs(order, groups, groups))
scatt_coeffs = ZERO
allocate(ScattDataTabular :: this % scatter)
else if (scatt_type == ANGLE_LEGENDRE) then
! Otherwise find the maximum scattering order
! Need to determine the maximum scattering order of all data in this material
mat_max_order = 0
do i = 1, mat % n_nuclides
if (nuclides(mat % nuclide(i)) % obj % order > mat_max_order) then
mat_max_order = nuclides(mat % nuclide(i)) % obj % order
end if
end do
! Now need to compare this material maximum scattering order with
! the problem wide max scatt order and use whichever is lower
order = min(mat_max_order, max_order)
this % order = order + 1
! Now we can allocate our scatt_coeffs object accordingly
allocate(scatt_coeffs(order + 1, groups, groups))
scatt_coeffs = ZERO
if (legendre_mu_points == 1) then
allocate(ScattDataLegendre :: this % scatter)
else
allocate(ScattDataTabular :: this % scatter)
end if
end if
! Allocate and initialize data within macro_xs(i_mat) object
allocate(this % total(groups))
this % total = ZERO
allocate(this % absorption(groups))
this % absorption = ZERO
if (get_fiss) then
allocate(this % fission(groups))
this % fission = ZERO
end if
if (get_kfiss) then
allocate(this % k_fission(groups))
this % k_fission = ZERO
end if
allocate(this % nu_fission(groups))
this % nu_fission = ZERO
allocate(this % chi(groups, groups))
this % chi = ZERO
allocate(temp_energy(groups, groups))
temp_energy = ZERO
allocate(temp_mult(groups, groups))
temp_mult = ZERO
allocate(this % scattxs(groups))
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
! Perform our operations which depend upon the type
select type(nuc => nuclides(mat % nuclide(i)) % obj)
type is (NuclideIso)
! Add contributions to total, absorption, and fission data (if necessary)
this % total = this % total + atom_density * nuc % total
this % absorption = this % absorption + &
atom_density * nuc % absorption
if (nuc % fissionable) then
if (allocated(nuc % chi)) then
do gin = 1, groups
do gout = 1, groups
this % chi(gout,gin) = this % chi(gout,gin) + atom_density * &
nuc % chi(gout) * nuc % nu_fission(gin,1)
end do
end do
this % nu_fission = this % nu_fission + atom_density * &
nuc % nu_fission(:,1)
else
this % chi = this % chi + atom_density * nuc % nu_fission
do gin = 1, groups
this % nu_fission(gin) = this % nu_fission(gin) + atom_density * &
sum(nuc % nu_fission(:,gin))
end do
end if
if (get_fiss) then
this % fission = this % fission + atom_density * nuc % fission
end if
if (get_kfiss) then
this % k_fission = this % k_fission + atom_density * nuc % k_fission
end if
end if
! Now time to do the scattering
do gin = 1, groups
do gout = 1, groups
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
! Transfer matrix
temp_energy(gout,gin) = temp_energy(gout,gin) + atom_density * &
sum(nuc % scatter(gout,gin,:))
! Determine the angular distribution
do imu = 1, order
scatt_coeffs(imu, gout, gin) = scatt_coeffs(imu, gout, gin) + &
nuc % scatter(gout,gin,imu) * &
atom_density
end do
else if (scatt_type == ANGLE_LEGENDRE) then
! Transfer matrix
temp_energy(gout,gin) = temp_energy(gout,gin) + atom_density * &
nuc % scatter(gout,gin,1)
! Determine the angular distribution coefficients so we can later
! expand do the complete distribution
do l = 1, min(nuc % order, order) + 1
scatt_coeffs(l, gout, gin) = scatt_coeffs(l, gout, gin) + &
nuc % scatter(gout,gin,l) * &
atom_density
end do
end if
! Multiplicity matrix
temp_mult(gout,gin) = temp_mult(gout,gin) + atom_density * &
nuc % mult(gout,gin)
end do
end do
type is (NuclideAngle)
error_code = 1
error_text = "Invalid Passing of NuclideAngle to MacroXSIso Object"
return
end select
end do
! Store the scattering xs
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
this % scattxs(:) = sum(sum(scatt_coeffs(:,:,:),dim=1),dim=1)
else if (scatt_type == ANGLE_LEGENDRE) then
this % scattxs(:) = sum(scatt_coeffs(1,:,:),dim=1)
end if
! Normalize the scatt_coeffs
do gin = 1, groups
do gout = 1, groups
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
norm = sum(scatt_coeffs(:,gout,gin))
else if (scatt_type == ANGLE_LEGENDRE) then
norm = scatt_coeffs(1,gout,gin)
end if
if (norm /= ZERO) then
scatt_coeffs(:, gout, gin) = scatt_coeffs(:, gout,gin) / norm
end if
end do
! Now normalize temp_energy (outgoing scattering energy probabilities)
norm = sum(temp_energy(:,gin))
if (norm > ZERO) then
temp_energy(:,gin) = temp_energy(:,gin) / norm
end if
end do
if (scatt_type == ANGLE_LEGENDRE .and. legendre_mu_points /= 1) then
call this % scatter % init(legendre_mu_points, temp_energy, temp_mult, &
scatt_coeffs)
else
call this % scatter % init(this % order, temp_energy, temp_mult, &
scatt_coeffs)
end if
! Now normalize chi
if (mat % fissionable) then
do gin = 1, groups
! Normalize Chi
norm = sum(this % chi(:,gin))
if (norm > ZERO) then
this % chi(:,gin) = this % chi(:,gin) / norm
end if
end do
end if
! Deallocate temporaries for the next material
deallocate(scatt_coeffs, temp_energy, temp_mult)
end subroutine macroxsiso_init
subroutine macroxsangle_init(this, mat, nuclides, groups, get_kfiss, get_fiss, &
max_order, scatt_type, legendre_mu_points, error_code, error_text)
class(MacroXSAngle), intent(inout) :: this ! The MacroXS to initialize
type(Material), pointer, intent(in) :: mat ! base material
type(NuclideMGContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
logical, intent(in) :: get_kfiss ! Should we get kfiss data?
logical, intent(in) :: get_fiss ! Should we get fiss data?
integer, intent(in) :: max_order ! Maximum requested order
integer, intent(in) :: scatt_type ! Legendre or Tabular Scatt?
integer, intent(in) :: legendre_mu_points ! Treat as Leg or Tabular?
integer, intent(inout) :: error_code ! Code signifying error
character(MAX_LINE_LEN), intent(inout) :: error_text ! Error message to print
integer :: i ! loop index over nuclides
integer :: gin, gout ! group indices
real(8) :: atom_density ! atom density of a nuclide
integer :: ipol, iazi, npol, nazi
integer :: imu
real(8) :: norm
integer :: mat_max_order, order, l
real(8), allocatable :: temp_mult(:,:,:,:)
real(8), allocatable :: temp_energy(:,:,:,:)
real(8), allocatable :: scatt_coeffs(:,:,:,:,:)
! Initialize error data
error_code = 0
error_text = ''
! Get the number of each polar and azi angles and make sure all the
! NuclideAngle types have the same number of these angles
npol = -1
nazi = -1
do i = 1, mat % n_nuclides
select type(nuc => nuclides(mat % nuclide(i)) % obj)
type is (NuclideAngle)
if (npol == -1) then
npol = nuc % n_pol
nazi = nuc % n_azi
allocate(this % polar(npol))
this % polar = nuc % polar
allocate(this % azimuthal(nazi))
this % azimuthal = nuc % azimuthal
else
if ((npol /= nuc % n_pol) .or. (nazi /= nuc % n_azi)) then
error_code = 1
error_text = "All Angular Data Must Be Same Length!"
end if
end if
end select
end do
! If we have tabular only data, then make sure all datasets have same size
if (scatt_type == ANGLE_HISTOGRAM) then
! Check all scattering data of same size
order = nuclides(mat % nuclide(1)) % obj % order
do i = 2, mat % n_nuclides
if (order /= nuclides(mat % nuclide(i)) % obj % order) then
error_code = 1
error_text = "All Histogram Scattering Entries Must Be Same Length!"
return
end if
end do
! Ok, got our order, store it
this % order = order
! Allocate stuff for later
allocate(scatt_coeffs(order, groups, groups, nazi, npol))
scatt_coeffs = ZERO
allocate(this % scatter(nazi, npol))
do ipol = 1, npol
do iazi = 1, nazi
allocate(ScattDataHistogram :: this % scatter(iazi, ipol) % obj)
end do
end do
else if (scatt_type == ANGLE_TABULAR) then
! Check all scattering data of same size
order = nuclides(mat % nuclide(1)) % obj % order
do i = 2, mat % n_nuclides
if (order /= nuclides(mat % nuclide(i)) % obj % order) then
error_code = 1
error_text = "All Tabular Scattering Entries Must Be Same Length!"
return
end if
end do
! Ok, got our order, store it
this % order = order
! Allocate stuff for later
allocate(scatt_coeffs(order, groups, groups, nazi, npol))
scatt_coeffs = ZERO
allocate(this % scatter(nazi, npol))
do ipol = 1, npol
do iazi = 1, nazi
allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj)
end do
end do
else if (scatt_type == ANGLE_LEGENDRE) then
! Otherwise find the maximum scattering order
! Need to determine the maximum scattering order of all data in this material
mat_max_order = 0
do i = 1, mat % n_nuclides
if (nuclides(mat % nuclide(i)) % obj % order > mat_max_order) then
mat_max_order = nuclides(mat % nuclide(i)) % obj % order
end if
end do
! Now need to compare this material maximum scattering order with
! the problem wide max scatt order and use whichever is lower
order = min(mat_max_order, max_order)
this % order = order + 1
! Now we can allocate our scatt_coeffs object accordingly
allocate(scatt_coeffs(order + 1, groups, groups, nazi, npol))
scatt_coeffs = ZERO
allocate(this % scatter(nazi, npol))
do ipol = 1, npol
do iazi = 1, nazi
if (legendre_mu_points == 1) then
allocate(ScattDataLegendre :: this % scatter(iazi, ipol) % obj)
else
allocate(ScattDataTabular :: this % scatter(iazi, ipol) % obj)
end if
end do
end do
end if
! Allocate and initialize data within macro_xs(i_mat) object
allocate(this % total(groups,nazi,npol))
this % total = ZERO
allocate(this % absorption(groups,nazi,npol))
this % absorption = ZERO
if (get_fiss) then
allocate(this % fission(groups,nazi,npol))
this % fission = ZERO
end if
if (get_kfiss) then
allocate(this % k_fission(groups,nazi,npol))
this % k_fission = ZERO
end if
allocate(this % nu_fission(groups,nazi,npol))
this % nu_fission = ZERO
allocate(this % chi(groups, groups, nazi, npol))
this % chi = ZERO
allocate(temp_energy(groups,groups,nazi,npol))
temp_energy = ZERO
allocate(temp_mult(groups,groups,nazi,npol))
temp_mult = ZERO
allocate(this % scattxs(groups,nazi,npol))
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! Copy atom density of nuclide in material
atom_density = mat % atom_density(i)
! Perform our operations which depend upon the type
select type(nuc => nuclides(mat % nuclide(i)) % obj)
type is (NuclideIso)
error_code = 1
error_text = "Invalid Passing of NuclideIso to MacroXSAngle Object"
return
type is (NuclideAngle)
! Add contributions to total, absorption, and fission data (if necessary)
this % total = this % total + atom_density * nuc % total
this % absorption = this % absorption + &
atom_density * nuc % absorption
if (nuc % fissionable) then
if (allocated(nuc % chi)) then
do gin = 1, groups
do gout = 1, groups
this % chi(gout,gin,:,:) = this % chi(gout,gin,:,:) + atom_density * &
nuc % chi(gout,:,:) * nuc % nu_fission(gin,1,:,:)
end do
end do
this % nu_fission = this % nu_fission + atom_density * &
nuc % nu_fission(:,1,:,:)
else
this % chi = this % chi + atom_density * nuc % nu_fission
do gin = 1, groups
this % nu_fission(gin,:,:) = this % nu_fission(gin,:,:) + atom_density * &
sum(nuc % nu_fission(:,gin,:,:),dim=1)
end do
end if
if (get_fiss) then
this % fission = this % fission + atom_density * nuc % fission
end if
if (get_kfiss) then
this % k_fission = this % k_fission + atom_density * nuc % k_fission
end if
end if
! Now time to do the scattering
do gin = 1, groups
do gout = 1, groups
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
! Transfer matrix
temp_energy(gout,gin,:,:) = temp_energy(gout,gin,:,:) + atom_density * &
sum(nuc % scatter(gout,gin,:,:,:),dim=1)
! Determine the angular distribution
do imu = 1, order
scatt_coeffs(imu,gout,gin,:,:) = scatt_coeffs(imu,gout,gin,:,:) + &
nuc % scatter(gout,gin,imu,:,:) * &
atom_density
end do
else if (scatt_type == ANGLE_LEGENDRE) then
! Transfer matrix
temp_energy(gout,gin,:,:) = temp_energy(gout,gin,:,:) + atom_density * &
nuc % scatter(gout,gin,1,:,:)
! Determine the angular distribution coefficients so we can later
! expand do the complete distribution
do l = 1, min(nuc % order, order) + 1
scatt_coeffs(l, gout, gin,:,:) = scatt_coeffs(l, gout, gin,:,:) + &
nuc % scatter(gout,gin,l,:,:) * &
atom_density
end do
end if
! Multiplicity matrix
temp_mult(gout,gin,:,:) = temp_mult(gout,gin,:,:) + atom_density * &
nuc % mult(gout,gin,:,:)
end do
end do
end select
end do
! Store the scattering xs
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
this % scattxs(:,:,:) = sum(sum(scatt_coeffs(:,:,:,:,:),dim=1),dim=1)
else if (scatt_type == ANGLE_LEGENDRE) then
this % scattxs(:,:,:) = sum(scatt_coeffs(1,:,:,:,:),dim=1)
end if
! Normalize the scatt_coeffs
do ipol = 1, npol
do iazi = 1, nazi
do gin = 1, groups
do gout = 1, groups
if (scatt_type == ANGLE_HISTOGRAM .or. scatt_type == ANGLE_TABULAR) then
norm = sum(scatt_coeffs(:,gout,gin,iazi,ipol))
else if (scatt_type == ANGLE_LEGENDRE) then
norm = scatt_coeffs(1,gout,gin,iazi,ipol)
end if
if (norm /= ZERO) then
scatt_coeffs(:,gout,gin,iazi,ipol) = &
scatt_coeffs(:,gout,gin,iazi,ipol) / norm
end if
end do
! Now normalize temp_energy (outgoing scattering energy probabilities)
norm = sum(temp_energy(:,gin,iazi,ipol))
if (norm > ZERO) then
temp_energy(:,gin,iazi,ipol) = temp_energy(:,gin,iazi,ipol) / norm
end if
end do
if (scatt_type == ANGLE_LEGENDRE .and. legendre_mu_points /= 1) then
call this % scatter(iazi, ipol) % obj % init(legendre_mu_points, &
temp_energy(:,:,iazi,ipol), temp_mult(:,:,iazi,ipol), &
scatt_coeffs(:,:,:,iazi,ipol))
else
call this % scatter(iazi, ipol) % obj % init(this % order, &
temp_energy(:,:,iazi,ipol), temp_mult(:,:,iazi,ipol), &
scatt_coeffs(:,:,:,iazi,ipol))
end if
end do
end do
! Now go through and normalize chi
if (mat % fissionable) then
do ipol = 1, npol
do iazi = 1, nazi
do gin = 1, groups
! Normalize Chi
norm = sum(this % chi(:,gin,iazi,ipol))
if (norm > ZERO) then
this % chi(:,gin,iazi,ipol) = this % chi(:,gin,iazi,ipol) / norm
end if
end do
end do
end do
end if
! Deallocate temporaries for the next material
deallocate(scatt_coeffs, temp_energy, temp_mult)
end subroutine macroxsangle_init
!===============================================================================
! MACROXS_*_GET_XS returns the requested data type
!===============================================================================
function macroxsiso_get_xs(this, g, xstype, gout, uvw) result(xs)
class(MacroXSIso), intent(in) :: this ! The MacroXS to initialize
integer, intent(in) :: g ! Incoming Energy group
character(*) , intent(in) :: xstype ! Type of xs requested
integer, optional, intent(in) :: gout ! Outgoing Energy group
real(8), optional, intent(in) :: uvw(3) ! Requested Angle
real(8) :: xs ! Requested x/s
select case(xstype)
case('total')
xs = this % total(g)
case('absorption')
xs = this % absorption(g)
case('fission')
xs = this % fission(g)
case('k_fission')
xs = this % k_fission(g)
case('nu_fission')
xs = this % nu_fission(g)
case('scatter')
xs = this % scattxs(g)
case('mult')
if (present(gout)) then
xs = this % scatter % mult(gout,g)
else
xs = sum(this % scatter % mult(:,g))
end if
end select
end function macroxsiso_get_xs
function macroxsangle_get_xs(this, g, xstype, gout,uvw) result(xs)
class(MacroXSAngle), intent(in) :: this ! The MacroXS to initialize
integer, intent(in) :: g ! Incoming Energy group
character(*) , intent(in) :: xstype ! Type of xs requested
integer, optional, intent(in) :: gout ! Outgoing Energy group
real(8), optional, intent(in) :: uvw(3) ! Requested Angle
real(8) :: xs ! Requested x/s
integer :: iazi, ipol
if (present(uvw)) then
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
select case(xstype)
case('total')
xs = this % total(g,iazi,ipol)
case('absorption')
xs = this % absorption(g,iazi,ipol)
case('fission')
xs = this % fission(g,iazi,ipol)
case('k_fission')
xs = this % k_fission(g,iazi,ipol)
case('nu_fission')
xs = this % nu_fission(g,iazi,ipol)
case('scatter')
xs = this % scattxs(g,iazi,ipol)
case('mult')
if (present(gout)) then
xs = this % scatter(iazi,ipol) % obj % mult(gout,g)
else
xs = sum(this % scatter(iazi,ipol) % obj % mult(:,g))
end if
end select
end if
end function macroxsangle_get_xs
!===============================================================================
! MACROXS_*_SAMPLE_FISSION_ENERGY samples the outgoing energy from a fission
! event
!===============================================================================
function macroxsiso_sample_fission_energy(this, gin, uvw) result(gout)
class(MacroXSIso), intent(in) :: this ! Data to work with
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
xi = prn()
prob = ZERO
gout = 0
do while (prob < xi)
gout = gout + 1
prob = prob + this % chi(gout,gin)
end do
end function macroxsiso_sample_fission_energy
function macroxsangle_sample_fission_energy(this, gin, uvw) result(gout)
class(MacroXSAngle), intent(in) :: this ! Data to work with
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
integer :: iazi, ipol
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
xi = prn()
prob = ZERO
gout = 0
do while (prob < xi)
gout = gout + 1
prob = prob + this % chi(gout,gin,iazi,ipol)
end do
end function macroxsangle_sample_fission_energy
!===============================================================================
! MACROXS*_SAMPLE_SCATTER Selects outgoing energy and angle after a scatter
! event
!===============================================================================
subroutine macroxsiso_sample_scatter(this, uvw, gin, gout, mu, wgt)
class(MacroXSIso), intent(in) :: this
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
call this % scatter % sample(gin, gout, mu, wgt)
end subroutine macroxsiso_sample_scatter
subroutine macroxsangle_sample_scatter(this, uvw, gin, gout, mu, wgt)
class(MacroXSAngle), intent(in) :: this
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
integer :: iazi, ipol ! Angular indices
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
call this % scatter(iazi,ipol) % obj % sample(gin,gout,mu,wgt)
end subroutine macroxsangle_sample_scatter
!===============================================================================
! MACROXS*_CALCULATE_XS determines the multi-group macroscopic cross sections
! for the material the particle is currently traveling through.
!===============================================================================
subroutine macroxsiso_calculate_xs(this, gin, uvw, xs)
class(MacroXSIso), intent(in) :: this
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs ! Resultant MacroXS Data
xs % total = this % total(gin)
xs % elastic = this % scattxs(gin)
xs % absorption = this % absorption(gin)
xs % nu_fission = this % nu_fission(gin)
end subroutine macroxsiso_calculate_xs
subroutine macroxsangle_calculate_xs(this, gin, uvw, xs)
class(MacroXSAngle), intent(in) :: this
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs ! Resultant MacroXS Data
integer :: iazi, ipol
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
xs % total = this % total(gin, iazi, ipol)
xs % elastic = this % scattxs(gin, iazi, ipol)
xs % absorption = this % absorption(gin, iazi, ipol)
xs % nu_fission = this % nu_fission(gin, iazi, ipol)
end subroutine macroxsangle_calculate_xs
end module macroxs_header

70
src/math.F90
View file

@ -557,6 +557,49 @@ contains
end function calc_rn
!===============================================================================
! EXPAND_HARMONIC expands a given series of real spherical harmonics
!===============================================================================
pure function expand_harmonic(data, order, uvw) result(val)
real(8), intent(in) :: data(:)
integer, intent(in) :: order
real(8), intent(in) :: uvw(3)
real(8) :: val
integer :: l, lm_lo, lm_hi
val = data(1)
lm_lo = 2
lm_hi = 4
do l = 1, order - 1
val = val + sqrt(TWO * real(l,8) + ONE) * &
dot_product(calc_rn(l,uvw), data(lm_lo:lm_hi))
lm_lo = lm_hi + 1
lm_hi = lm_lo + 2 * (l + 1)
end do
end function expand_harmonic
!===============================================================================
! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
! and the value of x
!===============================================================================
pure function evaluate_legendre(data, x) result(val)
real(8), intent(in) :: data(:)
real(8), intent(in) :: x
real(8) :: val
integer :: l
val = HALF * data(1)
do l = 1, size(data) - 1
val = val + (real(l,8) + HALF) * data(l + 1) * calc_pn(l,x)
end do
end function evaluate_legendre
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
@ -659,4 +702,31 @@ contains
end function watt_spectrum
!===============================================================================
! find_angle finds the closest angle on the data grid and returns that index
!===============================================================================
pure subroutine find_angle(polar, azimuthal, uvw, i_azi, i_pol)
real(8), intent(in) :: polar(:) ! Polar angles [0,pi]
real(8), intent(in) :: azimuthal(:) ! Azi. angles [-pi,pi]
real(8), intent(in) :: uvw(3) ! Direction of motion
integer, intent(inout) :: i_pol ! Closest polar bin
integer, intent(inout) :: i_azi ! Closest azi bin
real(8) :: my_pol, my_azi, dangle
! Convert uvw to polar and azi
my_pol = acos(uvw(3))
my_azi = atan2(uvw(2), uvw(1))
! Search for equi-binned angles
dangle = PI / real(size(polar),8)
i_pol = floor(my_pol / dangle + ONE)
dangle = TWO * PI / real(size(azimuthal),8)
i_azi = floor((my_azi + PI) / dangle + ONE)
end subroutine find_angle
end module math

7
src/mesh.F90
View file

@ -3,7 +3,6 @@ module mesh
use constants
use global
use mesh_header
use particle_header, only: Particle
use search, only: binary_search
#ifdef MPI
@ -160,10 +159,10 @@ contains
sites_outside)
type(RegularMesh), pointer :: m ! mesh to count sites
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
! cell and energy group
real(8), intent(in), optional :: energies(:) ! energy grid to search
real(8), intent(in), optional :: energies(:) ! energy grid to search
integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc)
logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?

246
src/mgxs_data.F90 Normal file
View file

@ -0,0 +1,246 @@
module mgxs_data
use constants
use error, only: fatal_error
use global
use macroxs_header
use material_header, only: Material
use nuclide_header
use output, only: write_message
use set_header, only: SetChar
use string, only: to_lower
use xml_interface
implicit none
contains
!===============================================================================
! READ_XS reads all the cross sections for the problem and stores them in
! nuclides and sab_tables arrays
!===============================================================================
subroutine read_mgxs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_listing ! index in xs_listings array
integer :: i_nuclide ! index in nuclides
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of isotope, e.g. U-235.03c
integer :: representation ! Data representation
type(Material), pointer :: mat
type(SetChar) :: already_read
type(Node), pointer :: doc => null()
type(Node), pointer :: node_xsdata
type(NodeList), pointer :: node_xsdata_list => null()
logical :: file_exists
character(MAX_LINE_LEN) :: temp_str
logical :: get_kfiss, get_fiss
integer :: l
! Check if cross_sections.xml exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
&// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Loading Cross Section Data...", 5)
! Parse cross_sections.xml file
call open_xmldoc(doc, path_cross_sections)
! Get node list of all <xsdata>
call get_node_list(doc, "xsdata", node_xsdata_list)
n_listings = get_list_size(node_xsdata_list)
! allocate arrays for ACE table storage and cross section cache
allocate(nuclides_MG(n_nuclides_total))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! Find out if we need fission & kappa fission
! (i.e., are there any SCORE_FISSION or SCORE_KAPPA_FISSION tallies?)
get_kfiss = .false.
get_fiss = .false.
do i = 1, n_tallies
do l = 1, tallies(i) % n_score_bins
if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then
get_kfiss = .true.
end if
if (tallies(i) % score_bins(l) == SCORE_FISSION) then
get_fiss = .true.
end if
end do
if (get_kfiss .and. get_fiss) exit
end do
! ==========================================================================
! READ ALL ACE CROSS SECTION TABLES
! Loop over all files
MATERIAL_LOOP: do i = 1, n_materials
mat => materials(i)
NUCLIDE_LOOP: do j = 1, mat % n_nuclides
name = mat % names(j)
if (.not. already_read % contains(name)) then
i_listing = xs_listing_dict % get_key(to_lower(name))
i_nuclide = mat % nuclide(j)
name = xs_listings(i_listing) % name
alias = xs_listings(i_listing) % alias
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i_listing, node_xsdata)
call write_message("Loading " // trim(name) // " Data...", 5)
! First find out the data representation
if (check_for_node(node_xsdata, "representation")) then
call get_node_value(node_xsdata, "representation", temp_str)
temp_str = trim(to_lower(temp_str))
if (temp_str == 'isotropic' .or. temp_str == 'iso') then
representation = MGXS_ISOTROPIC
else if (temp_str == 'angle') then
representation = MGXS_ANGLE
else
call fatal_error("Invalid Data Representation!")
end if
else
! Default to isotropic representation
representation = MGXS_ISOTROPIC
end if
! Now allocate accordingly
select case(representation)
case(MGXS_ISOTROPIC)
allocate(NuclideIso :: nuclides_MG(i_nuclide) % obj)
case(MGXS_ANGLE)
allocate(NuclideAngle :: nuclides_MG(i_nuclide) % obj)
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % init(node_xsdata, energy_groups, &
get_kfiss, get_fiss)
! Keep track of what listing is associated with this nuclide
nuclides_MG(i_nuclide) % obj % listing = i_listing
! Add name and alias to dictionary
call already_read % add(name)
call already_read % add(alias)
end if
end do NUCLIDE_LOOP
end do MATERIAL_LOOP
! Avoid some valgrind leak errors
call already_read % clear()
! Loop around material
MATERIAL_LOOP3: do i = 1, n_materials
! Get material
mat => materials(i)
! Loop around nuclides in material
NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! Is this fissionable?
if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
mat % fissionable = .true.
end if
if (mat % fissionable) then
exit NUCLIDE_LOOP2
end if
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
end subroutine read_mgxs
!===============================================================================
! SAME_NUCLIDEMG_LIST creates a linked list for each nuclide containing the
! indices in the nuclides array of all other instances of that nuclide. For
! example, the same nuclide may exist at multiple temperatures resulting
! in multiple entries in the nuclides array for a single zaid number.
!===============================================================================
subroutine same_nuclidemg_list()
integer :: i ! index in nuclides array
integer :: j ! index in nuclides array
do i = 1, n_nuclides_total
do j = 1, n_nuclides_total
if (nuclides_MG(i) % obj % zaid == nuclides_MG(j) % obj % zaid) then
call nuclides_MG(i) % obj % nuc_list % push_back(j)
end if
end do
end do
end subroutine same_nuclidemg_list
!===============================================================================
! CREATE_MACRO_XS generates the macroscopic x/s from the microscopic input data
!===============================================================================
subroutine create_macro_xs()
integer :: i_mat ! index in materials array
integer :: i ! loop index over nuclides
integer :: l ! Loop over score bins
type(Material), pointer :: mat ! current material
logical :: get_kfiss, get_fiss
integer :: error_code
character(MAX_LINE_LEN) :: error_text
integer :: scatt_type
integer :: legendre_mu_points
! Find out if we need fission & kappa fission
! (i.e., are there any SCORE_FISSION or SCORE_KAPPA_FISSION tallies?)
get_kfiss = .false.
get_fiss = .false.
do i = 1, n_tallies
do l = 1, tallies(i) % n_score_bins
if (tallies(i) % score_bins(l) == SCORE_KAPPA_FISSION) then
get_kfiss = .true.
end if
if (tallies(i) % score_bins(l) == SCORE_FISSION) then
get_fiss = .true.
end if
end do
if (get_kfiss .and. get_fiss) &
exit
end do
allocate(macro_xs(n_materials))
do i_mat = 1, n_materials
mat => materials(i_mat)
! Check to see how our nuclides are represented
! Force all to be the same type
! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs)
! At the same time, we will find the scattering type, as that will dictate
! how we allocate the scatter object within macroxs
legendre_mu_points = nuclides_MG(mat % nuclide(1)) % obj % legendre_mu_points
scatt_type = nuclides_MG(mat % nuclide(1)) % obj % scatt_type
select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
type is (NuclideIso)
allocate(MacroXSIso :: macro_xs(i_mat) % obj)
type is (NuclideAngle)
allocate(MacroXSAngle :: macro_xs(i_mat) % obj)
end select
call macro_xs(i_mat) % obj % init(mat, nuclides_MG, energy_groups, &
get_kfiss, get_fiss, max_order, &
scatt_type, legendre_mu_points, &
error_code, error_text)
! Handle any errors
if (error_code /= 0) call fatal_error(trim(error_text))
end do
end subroutine create_macro_xs
end module mgxs_data

1128
src/nuclide_header.F90 Normal file
View file

File diff suppressed because it is too large Load diff

210
src/output.F90
View file

@ -2,7 +2,7 @@ module output
use, intrinsic :: ISO_FORTRAN_ENV
use ace_header, only: Nuclide, Reaction, UrrData
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning
@ -12,8 +12,10 @@ module output
use math, only: t_percentile
use mesh_header, only: RegularMesh
use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices
use nuclide_header
use particle_header, only: LocalCoord, Particle
use plot_header
use sab_header, only: SAlphaBeta
use string, only: to_upper, to_str
use tally_header, only: TallyObject
@ -100,9 +102,9 @@ contains
!===============================================================================
subroutine header(msg, unit, level)
character(*), intent(in) :: msg ! header message
integer, intent(in), optional :: unit ! unit to write to
integer, intent(in), optional :: level ! specified header level
character(*), intent(in) :: msg ! header message
integer, intent(in), optional :: unit ! unit to write to
integer, intent(in), optional :: level ! specified header level
integer :: n ! number of = signs on left
integer :: m ! number of = signs on right
@ -305,175 +307,16 @@ contains
! Display weight, energy, grid index, and interpolation factor
write(ou,*) ' Weight = ' // to_str(p % wgt)
write(ou,*) ' Energy = ' // to_str(p % E)
if (run_CE) then
write(ou,*) ' Energy = ' // to_str(p % E)
else
write(ou,*) ' Energy Group = ' // to_str(p % g)
end if
write(ou,*) ' Delayed Group = ' // to_str(p % delayed_group)
write(ou,*)
end subroutine print_particle
!===============================================================================
! PRINT_NUCLIDE displays information about a continuous-energy neutron
! cross_section table and its reactions and secondary angle/energy distributions
!===============================================================================
subroutine print_nuclide(nuc, unit)
type(Nuclide), intent(in) :: nuc
integer, intent(in), optional :: unit
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
integer :: size_xs ! memory used for cross-sections (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
type(UrrData), pointer :: urr
! 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(nuc % name)
write(unit_,*) ' zaid = ' // trim(to_str(nuc % zaid))
write(unit_,*) ' awr = ' // trim(to_str(nuc % awr))
write(unit_,*) ' kT = ' // trim(to_str(nuc % kT))
write(unit_,*) ' # of grid points = ' // trim(to_str(nuc % n_grid))
write(unit_,*) ' Fissionable = ', nuc % fissionable
write(unit_,*) ' # of fission reactions = ' // trim(to_str(nuc % n_fission))
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, nuc % n_reaction
associate (rxn => nuc % 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 + (nuc % n_grid - rxn%threshold + 1) * 8
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
! nu-fission)
size_xs = 8 * nuc % n_grid * 4
! Write information about URR probability tables
size_urr = 0
if (nuc % urr_present) then
urr => nuc % 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 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 print_nuclide
!===============================================================================
! 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 print_sab_table(sab, unit)
type(SAlphaBeta), intent(in) :: sab
integer, intent(in), optional :: unit
integer :: size_sab ! memory used by S(a,b) table
integer :: unit_ ! unit to write to
integer :: i ! Loop counter for parsing through sab % 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(sab % name)
write(unit_,'(A)',advance="no") ' zaids = '
! Initialize the counter based on the above string
char_count = 11
do i = 1, sab % 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 < sab % n_zaid) then
! Include a comma
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", "
char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2
else
! Don't include a comma, since we are all done
write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i)))
end if
end do
write(unit_,*) "" ! Move to next line
write(unit_,*) ' awr = ' // trim(to_str(sab % awr))
write(unit_,*) ' kT = ' // trim(to_str(sab % kT))
! Inelastic data
write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_in))
write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_out))
write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_mu))
write(unit_,*) ' Threshold for Inelastic = ' // &
trim(to_str(sab % threshold_inelastic))
! Elastic data
if (sab % n_elastic_e_in > 0) then
write(unit_,*) ' # of Incoming Energies (Elastic) = ' // &
trim(to_str(sab % n_elastic_e_in))
write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // &
trim(to_str(sab % n_elastic_mu))
write(unit_,*) ' Threshold for Elastic = ' // &
trim(to_str(sab % threshold_elastic))
end if
! Determine memory used by S(a,b) table and write out
size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * &
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes'
! Blank line at end
write(unit_,*)
end subroutine print_sab_table
!===============================================================================
! 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
@ -486,8 +329,6 @@ contains
integer :: i ! loop index
integer :: unit_xs ! cross_sections.out file unit
character(MAX_FILE_LEN) :: path ! path of summary file
type(Nuclide), pointer :: nuc
type(SAlphaBeta), pointer :: sab
! Create filename for log file
path = trim(path_output) // "cross_sections.out"
@ -498,21 +339,22 @@ contains
! Write header
call header("CROSS SECTION TABLES", unit=unit_xs)
NUCLIDE_LOOP: do i = 1, n_nuclides_total
! Get pointer to nuclide
nuc => nuclides(i)
if (run_CE) then
NUCLIDE_LOOP: do i = 1, n_nuclides_total
! Print information about nuclide
call nuclides(i) % print(unit=unit_xs)
end do NUCLIDE_LOOP
! Print information about nuclide
call print_nuclide(nuc, unit=unit_xs)
end do NUCLIDE_LOOP
SAB_TABLES_LOOP: do i = 1, n_sab_tables
! Get pointer to S(a,b) table
sab => sab_tables(i)
! Print information about S(a,b) table
call print_sab_table(sab, unit=unit_xs)
end do SAB_TABLES_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
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
end if
! Close cross section summary file
close(unit_xs)

35
src/particle_header.F90
View file

@ -47,9 +47,14 @@ module particle_header
integer :: n_coord ! number of current coordinates
type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels
! Energy Data
real(8) :: E ! post-collision energy
real(8) :: last_E ! pre-collision energy
integer :: g ! post-collision energy group (MG only)
integer :: last_g ! pre-collision energy group (MG only)
! Other physical data
real(8) :: wgt ! particle weight
real(8) :: E ! energy
real(8) :: mu ! angle of scatter
logical :: alive ! is particle alive?
@ -57,7 +62,6 @@ module particle_header
real(8) :: last_xyz(3) ! previous coordinates
real(8) :: last_uvw(3) ! previous direction coordinates
real(8) :: last_wgt ! pre-collision particle weight
real(8) :: last_E ! pre-collision energy
real(8) :: absorb_wgt ! weight absorbed for survival biasing
! What event last took place
@ -128,6 +132,7 @@ contains
this % fission = .false.
this % delayed_group = 0
this % n_delayed_bank(:) = 0
this % g = 1
! Set up base level coordinates
this % coord(1) % universe = BASE_UNIVERSE
@ -174,9 +179,11 @@ contains
! fission, or simply as a secondary particle.
!===============================================================================
subroutine initialize_from_source(this, src)
class(Particle), intent(inout) :: this
type(Bank), intent(in) :: src
subroutine initialize_from_source(this, src, run_CE, energy_bin_avg)
class(Particle), intent(inout) :: this
type(Bank), intent(in) :: src
logical, intent(in) :: run_CE
real(8), allocatable, intent(in) :: energy_bin_avg(:)
! set defaults
call this % initialize()
@ -188,8 +195,14 @@ contains
this % coord(1) % uvw = src % uvw
this % last_xyz = src % xyz
this % last_uvw = src % uvw
this % E = src % E
this % last_E = src % E
if (run_CE) then
this % E = src % E
else
this % g = int(src % E)
this % last_g = int(src % E)
this % E = energy_bin_avg(this % g)
end if
this % last_E = src % E
end subroutine initialize_from_source
@ -198,10 +211,11 @@ contains
! the secondary bank and increments the number of sites in the secondary bank.
!===============================================================================
subroutine create_secondary(this, uvw, type)
subroutine create_secondary(this, uvw, type, run_CE)
class(Particle), intent(inout) :: this
real(8), intent(in) :: uvw(3)
integer, intent(in) :: type
logical, intent(in) :: run_CE
integer :: n
@ -215,8 +229,11 @@ contains
this % secondary_bank(n) % wgt = this % wgt
this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz
this % secondary_bank(n) % uvw(:) = uvw
this % secondary_bank(n) % E = this % E
this % n_secondary = n
this % secondary_bank(this % n_secondary) % E = this % E
if (.not. run_CE) then
this % secondary_bank(this % n_secondary) % E = real(this % g, 8)
end if
end subroutine create_secondary

2
src/particle_restart.F90
View file

@ -97,6 +97,7 @@ contains
call read_dataset(file_id, 'id', p%id)
call read_dataset(file_id, 'weight', p%wgt)
call read_dataset(file_id, 'energy', p%E)
call read_dataset(file_id, 'energy_group', p%g)
call read_dataset(file_id, 'xyz', p%coord(1)%xyz)
call read_dataset(file_id, 'uvw', p%coord(1)%uvw)
@ -105,6 +106,7 @@ contains
p%last_xyz = p%coord(1)%xyz
p%last_uvw = p%coord(1)%uvw
p%last_E = p%E
p%last_g = p%g
! Close hdf5 file
call file_close(file_id)

6
src/particle_restart_write.F90
View file

@ -1,10 +1,10 @@
module particle_restart_write
use bank_header, only: Bank
use bank_header, only: Bank
use global
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
use particle_header, only: Particle
use string, only: to_str
use hdf5

64
src/physics.F90
View file

@ -1,6 +1,6 @@
module physics
use ace_header, only: Nuclide, Reaction
use ace_header, only: Reaction
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
@ -9,15 +9,17 @@ module physics
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use math, only: rotate_angle, maxwell_spectrum, watt_spectrum
use math
use mesh, only: get_mesh_indices
use nuclide_header
use output, only: write_message
use particle_header, only: Particle
use particle_restart_write, only: write_particle_restart
use physics_common
use random_lcg, only: prn
use search, only: binary_search
use string, only: to_str
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use string, only: to_str
implicit none
@ -73,7 +75,7 @@ contains
integer :: i_nuclide ! index in nuclides array
integer :: i_nuc_mat ! index in material's nuclides array
integer :: i_reaction ! index in nuc % reactions array
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat)
@ -196,7 +198,7 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -242,7 +244,6 @@ contains
!===============================================================================
subroutine absorption(p, i_nuclide)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
@ -278,33 +279,11 @@ contains
end subroutine absorption
!===============================================================================
! RUSSIAN_ROULETTE
!===============================================================================
subroutine russian_roulette(p)
type(Particle), intent(inout) :: p
if (p % wgt < weight_cutoff) then
if (prn() < p % wgt / weight_survive) then
p % wgt = weight_survive
p % last_wgt = p % wgt
else
p % wgt = ZERO
p % last_wgt = ZERO
p % alive = .false.
end if
end if
end subroutine russian_roulette
!===============================================================================
! SCATTER
!===============================================================================
subroutine scatter(p, i_nuclide, i_nuc_mat)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_nuc_mat
@ -317,7 +296,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
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
! copy incoming direction
uvw_old(:) = p % coord(1) % uvw
@ -432,7 +411,7 @@ contains
real(8) :: v_cm(3) ! velocity of center-of-mass
real(8) :: v_t(3) ! velocity of target nucleus
real(8) :: uvw_cm(3) ! directional cosines in center-of-mass
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
! get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -501,7 +480,6 @@ contains
!===============================================================================
subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu)
integer, intent(in) :: i_nuclide ! index in micro_xs
integer, intent(in) :: i_sab ! index in sab_tables
real(8), intent(inout) :: E ! incoming/outgoing energy
@ -759,7 +737,7 @@ contains
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T
type(NuclideCE), 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
@ -1004,10 +982,10 @@ contains
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
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)
type(NuclideCE), 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) :: kT ! equilibrium temperature of target in MeV
real(8) :: awr ! target/neutron mass ratio
@ -1092,7 +1070,7 @@ contains
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide), pointer :: nuc
type(NuclideCE), pointer :: nuc
! Get pointers
nuc => nuclides(i_nuclide)
@ -1202,8 +1180,8 @@ contains
function sample_fission_energy(nuc, rxn, p) result(E_out)
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing energy of fission neutron
@ -1315,9 +1293,9 @@ contains
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, p)
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
type(NuclideCE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
integer :: i ! loop index
real(8) :: E ! energy in lab (incoming/outgoing)
@ -1365,7 +1343,7 @@ contains
p % wgt = yield * p % wgt
else
do i = 1, rxn % multiplicity - 1
call p % create_secondary(p % coord(1) % uvw, NEUTRON)
call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.True.)
end do
end if

33
src/physics_common.F90 Normal file
View file

@ -0,0 +1,33 @@
module physics_common
use constants
use global, only: weight_cutoff, weight_survive
use particle_header, only: Particle
use random_lcg, only: prn
implicit none
contains
!===============================================================================
! RUSSIAN_ROULETTE
!===============================================================================
subroutine russian_roulette(p)
type(Particle), intent(inout) :: p
if (p % wgt < weight_cutoff) then
if (prn() < p % wgt / weight_survive) then
p % wgt = weight_survive
p % last_wgt = p % wgt
else
p % wgt = ZERO
p % last_wgt = ZERO
p % alive = .false.
end if
end if
end subroutine russian_roulette
end module physics_common

271
src/physics_mg.F90 Normal file
View file

@ -0,0 +1,271 @@
module physics_mg
! This module contains the multi-group specific physics routines so as to not
! hinder performance of the CE versions with multiple if-thens.
use constants
use error, only: fatal_error, warning
use global
use macroxs_header, only: MacroXS, MacroXSContainer
use material_header, only: Material
use math, only: rotate_angle
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: Particle
use particle_restart_write, only: write_particle_restart
use physics_common
use random_lcg, only: prn
use scattdata_header
use string, only: to_str
implicit none
contains
!===============================================================================
! COLLISION_MG samples a nuclide and reaction and then calls the appropriate
! routine for that reaction
!===============================================================================
subroutine collision_mg(p)
type(Particle), intent(inout) :: p
! Store pre-collision particle properties
p % last_wgt = p % wgt
p % last_g = p % g
p % last_E = p % E
p % last_uvw = p % coord(1) % uvw
! Add to collision counter for particle
p % n_collision = p % n_collision + 1
! Sample nuclide/reaction for the material the particle is in
call sample_reaction(p)
! Display information about collision
if (verbosity >= 10 .or. trace) then
call write_message(" " // "Energy Group = " // trim(to_str(p % g)))
end if
end subroutine collision_mg
!===============================================================================
! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
! each nuclide within a material and then samples a reaction for that nuclide
! and calls the appropriate routine to process the physics. Note that there is
! special logic when suvival biasing is turned on since fission and
! disappearance are treated implicitly.
!===============================================================================
subroutine sample_reaction(p)
type(Particle), intent(inout) :: p
type(Material), pointer :: mat
mat => materials(p % material)
! Create fission bank sites. Note that while a fission reaction is sampled,
! it never actually "happens", i.e. the weight of the particle does not
! change when sampling fission sites. The following block handles all
! absorption (including fission)
if (mat % fissionable) then
if (run_mode == MODE_EIGENVALUE) then
call create_fission_sites(p, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE) then
call create_fission_sites(p, p % secondary_bank, p % n_secondary)
end if
end if
! If survival biasing is being used, the following subroutine adjusts the
! weight of the particle. Otherwise, it checks to see if absorption occurs
if (material_xs % absorption > ZERO) then
call absorption(p)
else
p % absorb_wgt = ZERO
end if
if (.not. p % alive) return
! Sample a scattering reaction and determine the secondary energy of the
! exiting neutron
call scatter(p)
! Play russian roulette if survival biasing is turned on
if (survival_biasing) then
call russian_roulette(p)
if (.not. p % alive) return
end if
end subroutine sample_reaction
!===============================================================================
! ABSORPTION
!===============================================================================
subroutine absorption(p)
type(Particle), intent(inout) :: p
if (survival_biasing) then
! Determine weight absorbed in survival biasing
p % absorb_wgt = (p % wgt * &
material_xs % absorption / material_xs % total)
! Adjust weight of particle by probability of absorption
p % wgt = p % wgt - p % absorb_wgt
p % last_wgt = p % wgt
! Score implicit absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
material_xs % nu_fission / material_xs % absorption
else
! See if disappearance reaction happens
if (material_xs % absorption > prn() * material_xs % total) then
! Score absorption estimate of keff
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % wgt * &
material_xs % nu_fission / material_xs % absorption
p % alive = .false.
p % event = EVENT_ABSORB
end if
end if
end subroutine absorption
!===============================================================================
! SCATTER
!===============================================================================
subroutine scatter(p)
type(Particle), intent(inout) :: p
call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
p % last_g, p % g, &
p % mu, p % wgt)
! Update energy value for downstream compatability (in tallying)
p % E = energy_bin_avg(p % g)
! Convert change in angle (mu) to new direction
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
! Set event component
p % event = EVENT_SCATTER
end subroutine scatter
!===============================================================================
! CREATE_FISSION_SITES determines the average total, prompt, and delayed
! neutrons produced from fission and creates appropriate bank sites.
!===============================================================================
subroutine create_fission_sites(p, bank_array, size_bank)
type(Particle), intent(inout) :: p
type(Bank), intent(inout) :: bank_array(:)
integer(8), intent(inout) :: size_bank
integer :: i ! loop index
integer :: nu ! actual number of neutrons produced
integer :: ijk(3) ! indices in ufs mesh
real(8) :: nu_t ! total nu
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
class(MacroXS), pointer :: xs
! Get Pointers
xs => macro_xs(p % material) % obj
! TODO: Heat generation from fission
! If uniform fission source weighting is turned on, we increase of decrease
! the expected number of fission sites produced
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh)
if (.not. in_mesh) then
call write_particle_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = ufs_mesh % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
else
weight = ONE
end if
! Determine expected number of neutrons produced
nu_t = p % wgt / keff * weight * &
material_xs % nu_fission / material_xs % total
! Sample number of neutrons produced
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
else
nu = int(nu_t) + 1
end if
! Check for bank size getting hit. For fixed source calculations, this is a
! fatal error. For eigenvalue calculations, it just means that k-effective
! was too high for a single batch.
if (size_bank + nu > size(bank_array)) then
if (run_mode == MODE_FIXEDSOURCE) then
call fatal_error("Secondary particle bank size limit reached. If you &
&are running a subcritical multiplication problem, k-effective &
&may be too close to one.")
else
if (master) call warning("Maximum number of sites in fission bank &
&reached. This can result in irreproducible results using different &
&numbers of processes/threads.")
end if
end if
! Bank source neutrons
if (nu == 0 .or. size_bank == size(bank_array)) return
p % fission = .true. ! Fission neutrons will be banked
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
! Bank source neutrons by copying particle data
bank_array(i) % xyz = p % coord(1) % xyz
! Set weight of fission bank site
bank_array(i) % wgt = ONE/weight
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
bank_array(i) % uvw(1) = mu
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
bank_array(i) % E = &
real(xs % sample_fission_energy(p % g, fission_bank(i) % uvw), 8)
end do
! increment number of bank sites
size_bank = min(size_bank + nu, int(size(bank_array),8))
! Store total weight banked for analog fission tallies
p % n_bank = nu
p % wgt_bank = nu/weight
end subroutine create_fission_sites
end module physics_mg

6
src/plot.F90
View file

@ -9,7 +9,7 @@ module plot
use mesh, only: get_mesh_indices
use mesh_header, only: RegularMesh
use output, only: write_message
use particle_header, only: Particle, LocalCoord
use particle_header, only: LocalCoord, Particle
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
@ -56,7 +56,7 @@ contains
subroutine position_rgb(p, pl, rgb, id)
type(Particle), intent(inout) :: p
type(Particle), intent(inout) :: p
type(ObjectPlot), pointer, intent(in) :: pl
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
@ -364,7 +364,7 @@ contains
real(8) :: ll(3) ! lower left starting point for each sweep direction
type(Particle) :: p
type(ProgressBar) :: progress
type(c_ptr) :: f_ptr
type(c_ptr) :: f_ptr
! compute voxel widths in each direction
vox = pl % width/dble(pl % pixels)

27
src/relaxng/materials.rnc
View file

@ -11,10 +11,10 @@ element materials {
} &
element nuclide {
(element name { xsd:string { maxLength = "7" } } |
(element name { xsd:string { maxLength = "7" } } |
attribute name { xsd:string { maxLength = "7" } }) &
(element xs { xsd:string { maxLength = "3" } } |
attribute xs { xsd:string { maxLength = "3" } })? &
(element xs { xsd:string { maxLength = "5" } } |
attribute xs { xsd:string { maxLength = "5" } })? &
(element scattering { ( "data" | "iso-in-lab" ) } |
attribute scattering { ( "data" | "iso-in-lab" ) })? &
(
@ -23,11 +23,18 @@ element materials {
)
}* &
element macroscopic {
(element name { xsd:string } |
attribute name { xsd:string }) &
(element xs { xsd:string { maxLength = "5" } } |
attribute xs { xsd:string { maxLength = "5" } })
}* &
element element {
(element name { xsd:string { maxLength = "2" } } |
(element name { xsd:string { maxLength = "2" } } |
attribute name { xsd:string { maxLength = "2" } }) &
(element xs { xsd:string { maxLength = "3" } } |
attribute xs { xsd:string { maxLength = "3" } })? &
(element xs { xsd:string { maxLength = "5" } } |
attribute xs { xsd:string { maxLength = "5" } })? &
(element scattering { ( "data" | "iso-in-lab" ) } |
attribute scattering { ( "data" | "iso-in-lab" ) })? &
(
@ -37,12 +44,12 @@ element materials {
}* &
element sab {
(element name { xsd:string { maxLength = "7" } } |
(element name { xsd:string { maxLength = "7" } } |
attribute name { xsd:string { maxLength = "7" } }) &
(element xs { xsd:string { maxLength = "3" } } |
attribute xs { xsd:string { maxLength = "3" } })?
(element xs { xsd:string { maxLength = "5" } } |
attribute xs { xsd:string { maxLength = "5" } })?
}*
}+ &
element default_xs { xsd:string { maxLength = "3" } }?
element default_xs { xsd:string { maxLength = "5" } }?
}

30
src/relaxng/materials.rng
View file

@ -118,6 +118,36 @@
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="macroscopic">
<interleave>
<choice>
<element name="name">
<data type="string">
<param name="maxLength">7</param>
</data>
</element>
<attribute name="name">
<data type="string">
<param name="maxLength">7</param>
</data>
</attribute>
</choice>
<choice>
<element name="xs">
<data type="string">
<param name="maxLength">3</param>
</data>
</element>
<attribute name="xs">
<data type="string">
<param name="maxLength">3</param>
</data>
</attribute>
</choice>
</interleave>
</element>
</zeroOrMore>
<zeroOrMore>
<element name="element">
<interleave>

61
src/relaxng/mg_cross_sections.rnc Normal file
View file

@ -0,0 +1,61 @@
element cross_sections {
element groups { xsd:int } &
element group_structure { list { xsd:double+ } } &
element inverse_velocities { list { xsd:double+ } }? &
element xsdata {
(element name { xsd:string { maxLength = "15" } } |
attribute name { xsd:string { maxLength = "15" } }) &
(element alias { xsd:string { maxLength = "15" } } |
attribute alias { xsd:string { maxLength = "15" } })? &
(element kT { xsd:double } | attribute kT { xsd:double })? &
(element fissionable { ( "true" | "false" ) } |
attribute fissionable { ( "true" | "false" ) }) &
(element representation { ( "isotropic" | "angle" ) } |
attribute representation { ( "isotropic" | "angle" ) })? &
(element num_azimuthal { xsd:positiveInteger } |
attribute num_azimuthal { xsd:positiveInteger })? &
(element num_polar { xsd:positiveInteger } |
attribute num_polar { xsd:positiveInteger })? &
(element scatt_type { ( "legendre" | "histogram" | "tabular" ) } |
attribute scatt_type { ( "legendre" | "histogram" | "tabular" ) })? &
(element order { xsd:positiveInteger } |
attribute order { xsd:positiveInteger }) &
element tabular_legendre {
(element enable { ( "true" | "false" ) } |
attribute enable { ( "true" | "false" ) })? &
(element num_points { xsd:positiveInteger } |
attribute num_points { xsd:positiveInteger })?
}? &
(element total { list { xsd:double+ } } |
attribute total { list { xsd:double+ } })? &
(element absorption { list { xsd:double+ } } |
attribute absorption { list { xsd:double+ } }) &
(element scatter { list { xsd:double+ } } |
attribute scatter { list { xsd:double+ } }) &
(element fission { list { xsd:double+ } } |
attribute fission { list { xsd:double+ } })? &
(element fission { list { xsd:double+ } } |
attribute fission { list { xsd:double+ } })? &
(element k_fission { list { xsd:double+ } } |
attribute k_fission { list { xsd:double+ } })? &
(element chi { list { xsd:double+ } } |
attribute chi { list { xsd:double+ } })? &
(element nu_fission { list { xsd:double+ } } |
attribute nu_fission { list { xsd:double+ } })?
}*
}

314
src/relaxng/mg_cross_sections.rng Normal file
View file

@ -0,0 +1,314 @@
<?xml version="1.0" encoding="UTF-8"?>
<element name="cross_sections" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<interleave>
<element name="groups">
<data type="int"/>
</element>
<element name="group_structure">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<optional>
<element name="inverse_velocities">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
</optional>
<zeroOrMore>
<element name="xsdata">
<interleave>
<choice>
<element name="name">
<data type="string">
<param name="maxLength">15</param>
</data>
</element>
<attribute name="name">
<data type="string">
<param name="maxLength">15</param>
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="alias">
<data type="string">
<param name="maxLength">15</param>
</data>
</element>
<attribute name="alias">
<data type="string">
<param name="maxLength">15</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="kT">
<data type="double"/>
</element>
<attribute name="kT">
<data type="double"/>
</attribute>
</choice>
</optional>
<choice>
<element name="fissionable">
<choice>
<value>true</value>
<value>false</value>
</choice>
</element>
<attribute name="fissionable">
<choice>
<value>true</value>
<value>false</value>
</choice>
</attribute>
</choice>
<optional>
<choice>
<element name="representation">
<choice>
<value>isotropic</value>
<value>angle</value>
</choice>
</element>
<attribute name="representation">
<choice>
<value>isotropic</value>
<value>angle</value>
</choice>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="num_azimuthal">
<data type="positiveInteger"/>
</element>
<attribute name="num_azimuthal">
<data type="positiveInteger"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="num_polar">
<data type="positiveInteger"/>
</element>
<attribute name="num_polar">
<data type="positiveInteger"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="scatt_type">
<choice>
<value>legendre</value>
<value>histogram</value>
<value>tabular</value>
</choice>
</element>
<attribute name="scatt_type">
<choice>
<value>legendre</value>
<value>histogram</value>
<value>tabular</value>
</choice>
</attribute>
</choice>
</optional>
<choice>
<element name="order">
<data type="positiveInteger"/>
</element>
<attribute name="order">
<data type="positiveInteger"/>
</attribute>
</choice>
<optional>
<element name="tabular_legendre">
<interleave>
<optional>
<choice>
<element name="enable">
<choice>
<value>true</value>
<value>false</value>
</choice>
</element>
<attribute name="enable">
<choice>
<value>true</value>
<value>false</value>
</choice>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="num_points">
<data type="positiveInteger"/>
</element>
<attribute name="num_points">
<data type="positiveInteger"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<choice>
<element name="total">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="total">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<choice>
<element name="absorption">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="absorption">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="scatter">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="scatter">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<optional>
<choice>
<element name="fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="k_fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="k_fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="chi">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="chi">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="nu_fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="nu_fission">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</zeroOrMore>
</interleave>
</element>

4
src/relaxng/settings.rnc
View file

@ -34,6 +34,8 @@ element settings {
element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? &
element energy_mode { ( "continuous-energy" | "ce" | "CE" | "multi-group" | "mg" | "MG" ) }? &
element entropy {
(element dimension { list { xsd:int+ } } |
attribute dimension { list { xsd:int+ } })? &
@ -45,6 +47,8 @@ element settings {
element log_grid_bins { xsd:positiveInteger }? &
element max_order { xsd:nonNegativeInteger }? &
element natural_elements { xsd:string { maxLength = "20" } }? &
element no_reduce { xsd:boolean }? &

17
src/relaxng/settings.rng
View file

@ -133,6 +133,18 @@
</choice>
</element>
</optional>
<optional>
<element name="energy_mode">
<choice>
<value>continuous-energy</value>
<value>ce</value>
<value>CE</value>
<value>multi-group</value>
<value>mg</value>
<value>MG</value>
</choice>
</element>
</optional>
<optional>
<element name="entropy">
<interleave>
@ -194,6 +206,11 @@
<data type="positiveInteger"/>
</element>
</optional>
<optional>
<element name="max_order">
<data type="nonNegativeInteger"/>
</element>
</optional>
<optional>
<element name="natural_elements">
<data type="string">

150
src/sab_header.F90 Normal file
View file

@ -0,0 +1,150 @@
module sab_header
use, intrinsic :: ISO_FORTRAN_ENV
use constants
use string, only: to_str
implicit none
!===============================================================================
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
! thermal scattering collisions which utilize a continuous secondary energy
! representation.
!===============================================================================
type DistEnergySab
integer :: n_e_out
real(8), allocatable :: e_out(:)
real(8), allocatable :: e_out_pdf(:)
real(8), allocatable :: e_out_cdf(:)
real(8), allocatable :: mu(:,:)
end type DistEnergySab
!===============================================================================
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
! of light isotopes such as water, graphite, Be, etc
!===============================================================================
type SAlphaBeta
character(10) :: 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
! threshold for S(a,b) treatment (usually ~4 eV)
real(8) :: threshold_inelastic
real(8) :: threshold_elastic = ZERO
! Inelastic scattering data
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
real(8), allocatable :: inelastic_e_out(:,:)
real(8), allocatable :: inelastic_mu(:,:,:)
! The following is used only if secondary_mode is 3
! The different implementation is necessary because the continuous
! representation has a variable number of outgoing energy points for each
! incoming energy
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
! Elastic scattering data
integer :: elastic_mode ! elastic mode (discrete/exact)
integer :: n_elastic_e_in ! # of incoming E for elastic
integer :: n_elastic_mu ! # of outgoing angles for elastic
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
real(8), allocatable :: elastic_mu(:,:)
contains
procedure :: print => print_sab_table
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 print_sab_table(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 print_sab_table
end module sab_header

489
src/scattdata_header.F90 Normal file
View file

@ -0,0 +1,489 @@
module scattdata_header
use constants
use error, only: fatal_error
use math
use random_lcg, only: prn
use search, only: binary_search
implicit none
!===============================================================================
! SCATTDATA contains all the data to describe the scattering energy and
! angular distribution
!===============================================================================
type, abstract :: ScattData
! p0 matrix on its own for sampling energy
real(8), allocatable :: energy(:,:) ! (Gout x Gin)
real(8), allocatable :: mult(:,:) ! (Gout x Gin)
real(8), allocatable :: data(:,:,:) ! (Order/Nmu x Gout x Gin)
contains
procedure(init_), deferred :: init ! Initializes ScattData
procedure(calc_f_), deferred :: calc_f ! Calculates f, given mu
procedure(sample_), deferred :: sample ! sample the scatter event
end type ScattData
abstract interface
subroutine init_(this, order, energy, mult, coeffs)
import ScattData
class(ScattData), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
end subroutine init_
pure function calc_f_(this, gin, gout, mu) result(f)
import ScattData
class(ScattData), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
end function calc_f_
subroutine sample_(this, gin, gout, mu, wgt)
import ScattData
class(ScattData), intent(in) :: this ! Scattering Object to Use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
end subroutine sample_
end interface
type, extends(ScattData) :: ScattDataLegendre
! Maximal value for rejection sampling from rectangle
real(8), allocatable :: max_val(:,:)
contains
procedure :: init => scattdatalegendre_init
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
end type ScattDataLegendre
type, extends(ScattData) :: ScattDataHistogram
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
contains
procedure :: init => scattdatahistogram_init
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
end type ScattDataHistogram
type, extends(ScattData) :: ScattDataTabular
real(8), allocatable :: mu(:) ! Mu bins
real(8) :: dmu ! Mu spacing
real(8), allocatable :: fmu(:,:,:) ! PDF of f(mu)
contains
procedure :: init => scattdatatabular_init
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
end type ScattDataTabular
!===============================================================================
! SCATTDATACONTAINER allocatable array for storing ScattData Objects (for angle)
!===============================================================================
type ScattDataContainer
class(ScattData), allocatable :: obj
end type ScattDataContainer
contains
!===============================================================================
! SCATTDATA_INIT builds the scattdata object
!===============================================================================
subroutine scattdata_init(this, order, energy, mult)
class(ScattData), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
integer :: groups
groups = size(energy, dim=1)
allocate(this % energy(groups, groups))
this % energy = energy
allocate(this % mult(groups, groups))
this % mult = mult
allocate(this % data(order, groups, groups))
this % data = ZERO
end subroutine scattdata_init
subroutine scattdatalegendre_init(this, order, energy, mult, coeffs)
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
real(8) :: dmu, mu, f
integer :: imu, Nmu, gout, gin, groups
call scattdata_init(this, order, energy, mult)
this % data = coeffs
groups = size(this % energy,dim=1)
allocate(this % max_val(groups, groups))
this % max_val = ZERO
! Step through the polynomial with fixed number of points to identify
! the maximal value.
Nmu = 1001
dmu = TWO / real(Nmu,8)
do imu = 1, Nmu
! Update mu. Do first and last seperate to avoid float errors
if (imu == 1) then
mu = -ONE
else if (imu == Nmu) then
mu = ONE
end if
mu = -ONE + real(imu - 1,8) * dmu
do gin = 1, groups
do gout = 1, groups
! Calculate probability
f = this % calc_f(gin,gout,mu)
! If this is a new max, store it.
if (f > this % max_val(gout,gin)) this % max_val(gout,gin) = f
end do
end do
end do
! Finally, since we may not have caught the exact max, add 10% margin
this % max_val = this % max_val * 1.1_8
end subroutine scattdatalegendre_init
subroutine scattdatahistogram_init(this, order, energy, mult, coeffs)
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
integer :: imu, gin, gout, groups
real(8) :: norm
groups = size(energy,dim=1)
call scattdata_init(this, order, energy, mult)
allocate(this % mu(order))
this % dmu = TWO / real(order,8)
this % mu(1) = -ONE
do imu = 2, order
this % mu(imu) = -ONE + real(imu - 1,8) * this % dmu
end do
! Best to integrate this histogram so we can avoid rejection sampling
do gin = 1, groups
do gout = 1, groups
if (energy(gout,gin) > ZERO) then
! Integrate the histogram
this % data(1,gout,gin) = this % dmu * coeffs(1,gout,gin)
do imu = 2, order
this % data(imu,gout,gin) = this % dmu * coeffs(imu,gout,gin) + &
this % data(imu-1,gout,gin)
end do
! Now make sure integral norms to zero
norm = this % data(order,gout,gin)
if (norm > ZERO) then
this % data(:,gout,gin) = this % data(:,gout,gin) / norm
end if
end if
end do
end do
end subroutine scattdatahistogram_init
subroutine scattdatatabular_init(this, order, energy, mult, coeffs)
class(ScattDataTabular), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
integer :: imu, gin, gout, groups
real(8) :: norm
logical :: legendre_flag
integer :: this_order
if (order < 0) then
legendre_flag = .true.
this_order = -1 * order
else
legendre_flag = .false.
this_order = order
end if
groups = size(energy,dim=1)
call scattdata_init(this, this_order, energy, mult)
allocate(this % mu(this_order))
this % dmu = TWO / real(this_order - 1)
do imu = 1, this_order - 1
this % mu(imu) = -ONE + real(imu - 1) * this % dmu
end do
this % mu(this_order) = ONE
! Best to integrate this histogram so we can avoid rejection sampling
allocate(this % fmu(this_order,groups,groups))
do gin = 1, groups
do gout = 1, groups
if (energy(gout,gin) > ZERO) then
if (legendre_flag) then
! Coeffs are legendre coeffs. Need to build f(mu) then integrate
! and store the integral in this % data
! Ensure the coeffs are normalized
norm = ONE / coeffs(1,gout,gin)
do imu = 1, this_order
this % fmu(imu,gout,gin) = evaluate_legendre(norm * coeffs(:,gout,gin), this % mu(imu))
! Force positivity
if (this % fmu(imu,gout,gin) < ZERO) then
this % fmu(imu,gout,gin) = ZERO
end if
end do
else
! Coeffs contain f(mu), put in f(mu) to save duplicate.
this % fmu(:,gout,gin) = this % data(:,gout,gin)
end if
! Re-normalize fmu for numerical integration issues and in case
! the negative fix-up introduced un-normalized data
norm = ZERO
do imu = 2, this_order
norm = norm + HALF * this % dmu * (this % fmu(imu-1,gout,gin) + this % fmu(imu,gout,gin))
end do
if (norm > ZERO) then
this % fmu(:,gout,gin) = this % fmu(:,gout,gin) / norm
end if
! Now create CDF from fmu with trapezoidal rule
this % data(1,gout,gin) = ZERO
do imu = 2, this_order - 1
this % data(imu,gout,gin) = this % data(imu-1,gout,gin) + &
HALF * this % dmu * (this % fmu(imu-1,gout,gin) + this % fmu(imu,gout,gin))
end do
this % data(this_order,gout,gin) = ONE
end if
end do
end do
end subroutine scattdatatabular_init
!===============================================================================
! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
!===============================================================================
pure function scattdatalegendre_calc_f(this, gin, gout, mu) result(f)
class(ScattDataLegendre), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
! Plug mu in to the legendre expansion and go from there
f = evaluate_legendre(this % data(:, gout, gin), mu)
end function scattdatalegendre_calc_f
pure function scattdatahistogram_calc_f(this, gin, gout, mu) result(f)
class(ScattDataHistogram), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
integer :: imu
! Find mu bin
imu = floor((mu + ONE)/ this % dmu + ONE)
! Adjust so interpolation works on the last bin if necessary
if (imu == size(this % data, dim=1)) then
imu = imu - 1
end if
! Use histogram interpolation to find f(mu)
f = this % data(imu, gout, gin)
end function scattdatahistogram_calc_f
pure function scattdatatabular_calc_f(this, gin, gout, mu) result(f)
class(ScattDataTabular), intent(in) :: this ! The ScattData to evaluate
integer, intent(in) :: gin ! Incoming Energy Group
integer, intent(in) :: gout ! Outgoing Energy Group
real(8), intent(in) :: mu ! Angle of interest
real(8) :: f ! Return value of f(mu)
integer :: imu
real(8) :: r
! Find mu bin
imu = floor((mu + ONE)/ this % dmu + ONE)
! Adjust so interpolation works on the last bin if necessary
if (imu == size(this % data, dim=1)) then
imu = imu - 1
end if
! ! Now interpolate to find f(mu)
r = (mu - this % mu(imu)) / (this % mu(imu + 1) - this % mu(imu))
f = (ONE - r) * this % data(imu, gout, gin) + &
r * this % data(imu + 1, gout, gin)
end function scattdatatabular_calc_f
!===============================================================================
! SCATTDATA*_SCATTER Samples the outgoing energy and change in angle.
!===============================================================================
subroutine scattdatalegendre_sample(this, gin, gout, mu, wgt)
class(ScattDataLegendre), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
real(8) :: u, f, M
integer :: samples
xi = prn()
prob = ZERO
gout = 0
do while (prob < xi)
gout = gout + 1
prob = prob + this % energy(gout,gin)
end do
! Now we can sample mu using the legendre representation of the thisering
! kernel in data(1:this % order)
! Do with rejection sampling
! Set maximal value
M = this % max_val(gout,gin)
samples = 0
do
mu = TWO * prn() - ONE
f = this % calc_f(gin,gout,mu)
if (f > ZERO) then
u = prn() * M
if (u <= f) then
exit
end if
end if
samples = samples + 1
if (samples > MAX_SAMPLE) then
call fatal_error("Maximum number of Legendre expansion samples reached!")
end if
end do
wgt = wgt * this % mult(gout,gin)
end subroutine scattdatalegendre_sample
subroutine scattdatahistogram_sample(this, gin, gout, mu, wgt)
class(ScattDataHistogram), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
integer :: imu
xi = prn()
prob = ZERO
gout = 0
do while (prob < xi)
gout = gout + 1
prob = prob + this % energy(gout,gin)
end do
xi = prn()
if (xi < this % data(1,gout,gin)) then
imu = 1
else
imu = binary_search(this % data(:,gout,gin), &
size(this % data(:,gout,gin)), xi)
end if
! Randomly select a mu in this bin.
mu = prn() * this % dmu + this % mu(imu)
wgt = wgt * this % mult(gout,gin)
end subroutine scattdatahistogram_sample
subroutine scattdatatabular_sample(this, gin, gout, mu, wgt)
class(ScattDataTabular), intent(in) :: this ! Scattering object to use
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
real(8) :: xi ! Our random number
real(8) :: prob ! Running probability
real(8) :: mu0, frac, mu1
real(8) :: c_k, c_k1, p0, p1
integer :: k, NP
xi = prn()
prob = ZERO
gout = 0
do while (prob < xi)
gout = gout + 1
prob = prob + this % energy(gout,gin)
end do
! determine outgoing cosine bin
NP = size(this % data(:,gout,gin))
xi = prn()
c_k = this % data(1,gout,gin)
do k = 1, NP - 1
c_k1 = this % data(k+1,gout,gin)
if (xi < c_k1) exit
c_k = c_k1
end do
! check to make sure k is <= NP - 1
k = min(k, NP - 1)
p0 = this % fmu(k,gout,gin)
mu0 = this % mu(k)
! Linear-linear interpolation to find mu value w/in bin.
p1 = this % fmu(k+1,gout,gin)
mu1 = this % mu(k+1)
frac = (p1 - p0)/(mu1 - mu0)
if (frac == ZERO) then
mu = mu0 + (xi - c_k)/p0
else
mu = mu0 + (sqrt(max(ZERO, p0*p0 + TWO*frac*(xi - c_k))) - p0)/frac
end if
if (mu <= -ONE) then
mu = -ONE
else if (mu >= ONE) then
mu = ONE
end if
wgt = wgt * this % mult(gout,gin)
end subroutine scattdatatabular_sample
end module scattdata_header

35
src/simulation.F90
View file

@ -4,26 +4,26 @@ module simulation
use mpi
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
use eigenvalue, only: join_bank_from_threads
#endif
use global
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: initialize_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
use random_lcg, only: set_particle_seed
use source, only: initialize_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
@ -131,7 +131,8 @@ contains
integer :: i
! set defaults
call p % initialize_from_source(source_bank(index_source))
call p % initialize_from_source(source_bank(index_source), run_CE, &
energy_bin_avg)
! set identifier for particle
p % id = work_index(rank) + index_source

53
src/source.F90
View file

@ -9,12 +9,13 @@ module source
use geometry_header, only: BASE_UNIVERSE
use global
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use math, only: maxwell_spectrum, watt_spectrum
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream
use state_point, only: read_source_bank, write_source_bank
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
@ -109,7 +110,7 @@ contains
integer, save :: num_resamples = 0 ! Number of resamples encountered
! Set weight to one by default
site%wgt = ONE
site % wgt = ONE
! Set the random number generator to the source stream.
call prn_set_stream(STREAM_SOURCE)
@ -117,10 +118,10 @@ contains
! Sample from among multiple source distributions
n_source = size(external_source)
if (n_source > 1) then
r(1) = prn()*sum(external_source(:)%strength)
r(1) = prn()*sum(external_source(:) % strength)
c = ZERO
do i = 1, n_source
c = c + external_source(i)%strength
c = c + external_source(i) % strength
if (r(1) < c) exit
end do
else
@ -131,14 +132,14 @@ contains
found = .false.
do while (.not.found)
! Set particle defaults
call p%initialize()
call p % initialize()
! Sample spatial distribution
site%xyz(:) = external_source(i)%space%sample()
site % xyz(:) = external_source(i) % space % sample()
! Fill p with needed data
p%coord(1)%xyz(:) = site%xyz
p%coord(1)%uvw(:) = [ ONE, ZERO, ZERO ]
p % coord(1) % xyz(:) = site % xyz
p % coord(1) % uvw(:) = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
@ -151,27 +152,27 @@ contains
end if
! Check if spatial site is in fissionable material
select type (space => external_source(i)%space)
select type (space => external_source(i) % space)
type is (SpatialBox)
if (space%only_fissionable) then
if (p%material == MATERIAL_VOID) then
if (space % only_fissionable) then
if (p % material == MATERIAL_VOID) then
found = .false.
elseif (.not. materials(p%material)%fissionable) then
elseif (.not. materials(p % material) % fissionable) then
found = .false.
end if
end if
end select
end do
call p%clear()
call p % clear()
! Sample angle
site%uvw(:) = external_source(i)%angle%sample()
site % uvw(:) = external_source(i) % angle % sample()
! Check for monoenergetic source above maximum neutron energy
select type (energy => external_source(i)%energy)
select type (energy => external_source(i) % energy)
type is (Discrete)
if (any(energy%x >= energy_max_neutron)) then
if (any(energy % x >= energy_max_neutron)) then
call fatal_error("Source energy above range of energies of at least &
&one cross section table")
end if
@ -179,14 +180,26 @@ contains
do
! Sample energy spectrum
site%E = external_source(i)%energy%sample()
site % E = external_source(i) % energy % sample()
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
if (site % E < energy_max_neutron) exit
end do
! Set delayed group
site%delayed_group = 0
site % delayed_group = 0
! If running in MG, convert site % E to group
if (.not. run_CE) then
if (site % E <= energy_bins(1)) then
site % E = real(1, 8)
else if (site % E > energy_bins(energy_groups + 1)) then
site % E = real(energy_groups, 8)
else
site % E = real(binary_search(energy_bins, energy_groups + 1, &
site % E), 8)
end if
end if
! Set the random number generator back to the tracking stream.
call prn_set_stream(STREAM_TRACKING)

19
src/state_point.F90
View file

@ -18,7 +18,7 @@ module state_point
use global
use hdf5_interface
use output, only: write_message, time_stamp
use string, only: to_str, zero_padded, count_digits
use string, only: to_str, count_digits, zero_padded
use tally_header, only: TallyObject
use mesh_header, only: RegularMesh
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
@ -93,6 +93,11 @@ contains
call write_dataset(file_id, "seed", seed)
! Write run information
if (run_CE) then
call write_dataset(file_id, "run_CE", 1)
else
call write_dataset(file_id, "run_CE", 0)
end if
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, "run_mode", "fixed source")
@ -647,6 +652,7 @@ contains
integer(HID_T) :: tally_group
real(8) :: real_array(3)
logical :: source_present
integer :: sp_run_CE
character(MAX_WORD_LEN) :: word
type(TallyObject), pointer :: tally
@ -674,6 +680,17 @@ contains
! Read and overwrite random number seed
call read_dataset(file_id, "seed", seed)
! It is not impossible for a state point to be generated from a CE run but
! to be loaded in to an MG run (or vice versa), check to prevent that.
call read_dataset(file_id, "run_CE", sp_run_CE)
if (sp_run_CE == 0 .and. run_CE) then
call fatal_error("State point file is from multi-group run but &
& current run is continous-energy!")
else if (sp_run_CE == 1 .and. .not. run_CE) then
call fatal_error("State point file is from continuous-energy run but &
& current run is multi-group!")
end if
! Read and overwrite run information except number of batches
call read_dataset(file_id, "run_mode", word)
select case(word)

15
src/summary.F90
View file

@ -1,6 +1,6 @@
module summary
use ace_header, only: Reaction, UrrData, Nuclide
use ace_header, only: Reaction, UrrData
use constants
use endf, only: reaction_name
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
@ -9,6 +9,7 @@ module summary
use hdf5_interface
use material_header, only: Material
use mesh_header, only: RegularMesh
use nuclide_header
use output, only: time_stamp
use surface_header
use string, only: to_str
@ -37,6 +38,12 @@ contains
! Write header information
call write_header(file_id)
if (run_CE) then
call write_dataset(file_id, "run_CE", 1)
else
call write_dataset(file_id, "run_CE", 0)
end if
! Write number of particles
call write_dataset(file_id, "n_particles", n_particles)
call write_dataset(file_id, "n_batches", n_batches)
@ -480,7 +487,11 @@ contains
! Copy ZAID for each nuclide to temporary array
allocate(nucnames(m%n_nuclides))
do j = 1, m%n_nuclides
i_list = nuclides(m%nuclide(j))%listing
if (run_CE) then
i_list = nuclides(m%nuclide(j))%listing
else
i_list = nuclides_MG(m%nuclide(j))%obj%listing
end if
nucnames(j) = xs_listings(i_list)%alias
end do

997
src/tally.F90
View file

File diff suppressed because it is too large Load diff

4
src/tally_header.F90
View file

@ -124,6 +124,10 @@ module tally_header
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Multi-Group Specific Information To Enable Rapid Tallying
logical :: energy_matches_groups = .false.
logical :: energyout_matches_groups = .false.
end type TallyObject
end module tally_header

4
src/track_output.F90
View file

@ -7,8 +7,8 @@ module track_output
use global
use hdf5_interface
use particle_header, only: Particle
use string, only: to_str
use particle_header, only: Particle
use string, only: to_str
use hdf5

75
src/tracking.F90
View file

@ -1,21 +1,23 @@
module tracking
use constants, only: MODE_EIGENVALUE
use cross_section, only: calculate_xs
use error, only: fatal_error, warning
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, check_cell_overlap
use geometry_header, only: Universe, BASE_UNIVERSE
use constants, only: MODE_EIGENVALUE
use cross_section, only: calculate_xs
use error, only: fatal_error, warning
use geometry, only: find_cell, distance_to_boundary, cross_surface, &
cross_lattice, check_cell_overlap
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use output, only: write_message
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use random_lcg, only: prn
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current
use track_output, only: initialize_particle_track, write_particle_track, &
add_particle_track, finalize_particle_track
use macroxs_header, only: MacroXS
use output, only: write_message
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use physics_mg, only: collision_mg
use random_lcg, only: prn
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current
use track_output, only: initialize_particle_track, write_particle_track, &
add_particle_track, finalize_particle_track
implicit none
@ -53,7 +55,9 @@ contains
total_weight = total_weight + p % wgt
! Force calculation of cross-sections by setting last energy to zero
micro_xs % last_E = ZERO
if (run_CE) then
micro_xs % last_E = ZERO
end if
! Prepare to write out particle track.
if (p % write_track) then
@ -79,11 +83,24 @@ contains
if (check_overlaps) call check_cell_overlap(p)
! Calculate microscopic and macroscopic cross sections -- note: 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)
! 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)
else
! Since the MGXS can be angle dependent, this needs to be done
! After every collision for the MGXS mode
if (p % material /= MATERIAL_VOID) then
call macro_xs(p % material) % obj % calculate_xs(p % g, &
p % coord(p % n_coord) % uvw, material_xs)
else
material_xs % total = ZERO
material_xs % elastic = ZERO
material_xs % absorption = ZERO
material_xs % nu_fission = ZERO
end if
end if
! Find the distance to the nearest boundary
call distance_to_boundary(p, d_boundary, surface_crossed, &
@ -105,8 +122,10 @@ contains
end do
! Score track-length tallies
if (active_tracklength_tallies % size() > 0) &
call score_tracklength_tally(p, distance)
if (active_tracklength_tallies % size() > 0) then
call score_tracklength_tally(p, distance)
end if
! Score track-length estimate of k-eff
if (run_mode == MODE_EIGENVALUE) then
@ -151,12 +170,15 @@ contains
! Clear surface component
p % surface = NONE
call collision(p)
if (run_CE) then
call collision(p)
else
call collision_mg(p)
end if
! Score collision estimator tallies -- this is done after a collision
! has occurred rather than before because we need information on the
! outgoing energy for any tallies with an outgoing energy filter
if (active_collision_tallies % size() > 0) call score_collision_tally(p)
if (active_analog_tallies % size() > 0) call score_analog_tally(p)
@ -200,7 +222,8 @@ contains
! Check for secondary particles if this particle is dead
if (.not. p % alive) then
if (p % n_secondary > 0) then
call p % initialize_from_source(p % secondary_bank(p % n_secondary))
call p % initialize_from_source(p % secondary_bank(p % n_secondary), &
run_CE, energy_bin_avg)
p % n_secondary = p % n_secondary - 1
n_event = 0

12
src/trigger.F90
View file

@ -6,12 +6,12 @@ module trigger
use constants
use global
use string, only: to_str
use output, only: warning, write_message
use mesh, only: mesh_indices_to_bin
use mesh_header, only: RegularMesh
use trigger_header, only: TriggerObject
use tally, only: TallyObject
use string, only: to_str
use output, only: warning, write_message
use mesh, only: mesh_indices_to_bin
use mesh_header, only: RegularMesh
use trigger_header, only: TriggerObject
use tally, only: TallyObject
implicit none

9859
tests/1d_mgxs.xml Normal file
View file

File diff suppressed because it is too large Load diff

89
tests/input_set.py
View file

@ -572,3 +572,92 @@ class InputSet(object):
plot.color = 'mat'
self.plots.add_plot(plot)
class MGInputSet(InputSet):
def build_default_materials_and_geometry(self):
# Define materials needed for 1D/1G slab problem
uo2_data = openmc.Macroscopic('uo2_iso', '71c')
uo2 = openmc.Material(name='UO2', material_id=1)
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
clad_data = openmc.Macroscopic('clad_ang_mu', '71c')
clad = openmc.Material(name='Clad', material_id=2)
clad.set_density('macro', 1.0)
clad.add_macroscopic(clad_data)
water_data = openmc.Macroscopic('lwtr_iso_mu', '71c')
water = openmc.Material(name='LWTR', material_id=3)
water.set_density('macro', 1.0)
water.add_macroscopic(water_data)
# Define the materials file.
self.materials.default_xs = '71c'
self.materials.add_materials((uo2, clad, water))
# Define surfaces.
# Assembly/Problem Boundary
left = openmc.XPlane(x0=0.0, surface_id=200,
boundary_type='reflective')
right = openmc.XPlane(x0=10.0, surface_id=201,
boundary_type='reflective')
bottom = openmc.YPlane(y0=0.0, surface_id=300,
boundary_type='reflective')
top = openmc.YPlane(y0=10.0, surface_id=301,
boundary_type='reflective')
down = openmc.ZPlane(z0=0.0, surface_id=0,
boundary_type='reflective')
fuel_clad_intfc = openmc.ZPlane(z0=2.0, surface_id=1)
clad_lwtr_intfc = openmc.ZPlane(z0=2.4, surface_id=2)
up = openmc.ZPlane(z0=5.0, surface_id=3,
boundary_type='reflective')
# Define cells
c1 = openmc.Cell(cell_id=1)
c1.region = +left & -right & +bottom & -top & +down & -fuel_clad_intfc
c1.fill = uo2
c2 = openmc.Cell(cell_id=2)
c2.region = +left & -right & +bottom & -top & +fuel_clad_intfc & -clad_lwtr_intfc
c2.fill = clad
c3 = openmc.Cell(cell_id=3)
c3.region = +left & -right & +bottom & -top & +clad_lwtr_intfc & -up
c3.fill = water
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cells((c1,c2,c3))
# Define the geometry file.
geometry = openmc.Geometry()
geometry.root_universe = root
self.geometry.geometry = geometry
def build_default_settings(self):
self.settings.batches = 10
self.settings.inactive = 5
self.settings.particles = 100
self.settings.source = Source(space=Box([0.0, 0.0, 0.0],
[10.0, 10.0, 2.0]))
self.settings.energy_mode = "multi-group"
self.settings.cross_sections = "../1d_mgxs.xml"
def build_defualt_plots(self):
plot = openmc.Plot()
plot.filename = 'mat'
plot.origin = (5.0, 5.0, 2.5)
plot.width = (2.5, 2.5)
plot.basis = 'xz'
plot.pixels = (3000, 3000)
plot.color = 'mat'
self.plots.add_plot(plot)

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

@ -0,0 +1 @@
04b4a5099f0097bbe02983c67dea691d0d0d4ece7fb7c264b9b2c29955baa9e870b6fa999480da08ead1e5a0c078ae33ce1b0a5c8594ad465aedf9bf3933e104

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

@ -0,0 +1,2 @@
k-combined:
1.045320E+00 5.851680E-02

17
tests/test_mg_basic/test_mg_basic.py Normal file
View file

@ -0,0 +1,17 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class MGBasicTestHarness(PyAPITestHarness):
def _build_inputs(self):
super(MGBasicTestHarness, self)._build_inputs()
if __name__ == '__main__':
harness = MGBasicTestHarness('statepoint.10.*', False, mg=True)
harness.main()

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

@ -0,0 +1 @@
abe20c626d613e73ccb1a3f8468ad1b9aecca528afa9e8131a411d754eb86b8ab64a6fb1fdc9c0b8b8158ff7c82f548de5912041bf035aa5a2d4532cfe0c9510

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

@ -0,0 +1,2 @@
k-combined:
1.083030E+00 1.855038E-02

85
tests/test_mg_max_order/test_mg_max_order.py Normal file
View file

@ -0,0 +1,85 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
from input_set import MGInputSet
import openmc
class MGNuclideInputSet(MGInputSet):
def build_default_materials_and_geometry(self):
# Define materials needed for 1D/1G slab problem
uo2_data = openmc.Macroscopic('uo2_iso', '71c')
uo2 = openmc.Material(name='UO2', material_id=1)
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
clad_data = openmc.Macroscopic('clad_iso', '71c')
clad = openmc.Material(name='Clad', material_id=2)
clad.set_density('macro', 1.0)
clad.add_macroscopic(clad_data)
water_data = openmc.Macroscopic('lwtr_iso', '71c')
water = openmc.Material(name='LWTR', material_id=3)
water.set_density('macro', 1.0)
water.add_macroscopic(water_data)
# Define the materials file.
self.materials.default_xs = '71c'
self.materials.add_materials((uo2, clad, water))
# Define surfaces.
# Assembly/Problem Boundary
left = openmc.XPlane(x0=0.0, surface_id=200,
boundary_type='reflective')
right = openmc.XPlane(x0=10.0, surface_id=201,
boundary_type='reflective')
bottom = openmc.YPlane(y0=0.0, surface_id=300,
boundary_type='reflective')
top = openmc.YPlane(y0=10.0, surface_id=301,
boundary_type='reflective')
down = openmc.ZPlane(z0=0.0, surface_id=0,
boundary_type='reflective')
fuel_clad_intfc = openmc.ZPlane(z0=2.0, surface_id=1)
clad_lwtr_intfc = openmc.ZPlane(z0=2.4, surface_id=2)
up = openmc.ZPlane(z0=5.0, surface_id=3,
boundary_type='reflective')
# Define cells
c1 = openmc.Cell(cell_id=1)
c1.region = +left & -right & +bottom & -top & +down & -fuel_clad_intfc
c1.fill = uo2
c2 = openmc.Cell(cell_id=2)
c2.region = +left & -right & +bottom & -top & +fuel_clad_intfc & -clad_lwtr_intfc
c2.fill = clad
c3 = openmc.Cell(cell_id=3)
c3.region = +left & -right & +bottom & -top & +clad_lwtr_intfc & -up
c3.fill = water
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cells((c1,c2,c3))
# Define the geometry file.
geometry = openmc.Geometry()
geometry.root_universe = root
self.geometry.geometry = geometry
class MGMaxOrderTestHarness(PyAPITestHarness):
def __init__(self, statepoint_name, tallies_present, mg=False):
PyAPITestHarness.__init__(self, statepoint_name, tallies_present)
self._input_set = MGNuclideInputSet()
def _build_inputs(self):
super(MGMaxOrderTestHarness, self)._build_inputs()
# Set P1 scattering
self._input_set.settings.max_order = 1
if __name__ == '__main__':
harness = MGMaxOrderTestHarness('statepoint.10.*', False, mg=True)
harness.main()

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

@ -0,0 +1 @@
c9f9e7211bfb2af58130bedfd64592d093b7bfa424953eba433ecf08940595a96b8de7a892f12d1ab465cebd8e5dd784114c1b1299b534ed329df92752c9ed1f

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

@ -0,0 +1,2 @@
k-combined:
1.380785E-01 5.556526E-03

84
tests/test_mg_nuclide/test_mg_nuclide.py Normal file
View file

@ -0,0 +1,84 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
from input_set import MGInputSet
import openmc
class MGNuclideInputSet(MGInputSet):
def build_default_materials_and_geometry(self):
# Define materials needed for 1D/1G slab problem
# This time do using nuclide, not macroscopic
uo2 = openmc.Material(name='UO2', material_id=1)
uo2.set_density('g/cm3', 1.0)
uo2.add_nuclide("uo2_iso", 1.0)
clad = openmc.Material(name='Clad', material_id=2)
clad.set_density('g/cm3', 1.0)
clad.add_nuclide("clad_ang_mu", 1.0)
water_data = openmc.Nuclide('lwtr_iso_mu', '71c')
water = openmc.Material(name='LWTR', material_id=3)
water.set_density('g/cm3', 1.0)
water.add_nuclide("lwtr_iso_mu", 1.0)
# Define the materials file.
self.materials.default_xs = '71c'
self.materials.add_materials((uo2, clad, water))
# Define surfaces.
# Assembly/Problem Boundary
left = openmc.XPlane(x0=0.0, surface_id=200,
boundary_type='reflective')
right = openmc.XPlane(x0=10.0, surface_id=201,
boundary_type='reflective')
bottom = openmc.YPlane(y0=0.0, surface_id=300,
boundary_type='reflective')
top = openmc.YPlane(y0=10.0, surface_id=301,
boundary_type='reflective')
down = openmc.ZPlane(z0=0.0, surface_id=0,
boundary_type='reflective')
fuel_clad_intfc = openmc.ZPlane(z0=2.0, surface_id=1)
clad_lwtr_intfc = openmc.ZPlane(z0=2.4, surface_id=2)
up = openmc.ZPlane(z0=5.0, surface_id=3,
boundary_type='reflective')
# Define cells
c1 = openmc.Cell(cell_id=1)
c1.region = +left & -right & +bottom & -top & +down & -fuel_clad_intfc
c1.fill = uo2
c2 = openmc.Cell(cell_id=2)
c2.region = +left & -right & +bottom & -top & +fuel_clad_intfc & -clad_lwtr_intfc
c2.fill = clad
c3 = openmc.Cell(cell_id=3)
c3.region = +left & -right & +bottom & -top & +clad_lwtr_intfc & -up
c3.fill = water
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cells((c1,c2,c3))
# Define the geometry file.
geometry = openmc.Geometry()
geometry.root_universe = root
self.geometry.geometry = geometry
class MGNuclideTestHarness(PyAPITestHarness):
def __init__(self, statepoint_name, tallies_present, mg=False):
PyAPITestHarness.__init__(self, statepoint_name, tallies_present)
self._input_set = MGNuclideInputSet()
def _build_inputs(self):
super(MGNuclideTestHarness, self)._build_inputs()
if __name__ == '__main__':
harness = MGNuclideTestHarness('statepoint.10.*', False, mg=True)
harness.main()

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

@ -0,0 +1 @@
ca8490e0e4549fed727ddc75b6d92cfe5162e11b905218a0afaa3ce2ee0763e2ff38074de27aaa678818624f49c5823650475dfa8f66f502a98fc03145399c0d

2906
tests/test_mg_tallies/results_true.dat Normal file
View file

File diff suppressed because it is too large Load diff

59
tests/test_mg_tallies/test_mg_tallies.py Normal file
View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
import openmc
class MGTalliesTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [17, 17, 1]
mesh.lower_left = [0.0, 0.0, 0.0]
mesh.upper_right = [21.42, 21.42, 100.0]
# Instantiate some tally filters
energy_filter = openmc.Filter(type='energy',
bins=[0.0, 20.0])
energyout_filter = openmc.Filter(type='energyout',
bins=[0.0, 20.0])
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
mat_filter = openmc.Filter(type='material', bins=[1,2,3])
tally1 = openmc.Tally(tally_id=1)
tally1.add_filter(mesh_filter)
tally1.add_score('total')
tally1.add_score('absorption')
tally1.add_score('flux')
tally1.add_score('fission')
tally1.add_score('nu-fission')
tally2 = openmc.Tally(tally_id=2)
tally2.add_filter(mat_filter)
tally2.add_filter(energy_filter)
tally2.add_filter(energyout_filter)
tally2.add_score('scatter')
tally2.add_score('nu-scatter')
self._input_set.tallies = openmc.TalliesFile()
self._input_set.tallies.add_mesh(mesh)
self._input_set.tallies.add_tally(tally1)
self._input_set.tallies.add_tally(tally2)
super(MGTalliesTestHarness, self)._build_inputs()
def _cleanup(self):
super(MGTalliesTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f): os.remove(f)
if __name__ == '__main__':
harness = MGTalliesTestHarness('statepoint.10.*', True, mg=True)
harness.main()

10
tests/testing_harness.py
View file

@ -12,7 +12,7 @@ import sys
import numpy as np
sys.path.insert(0, os.path.join(os.pardir, os.pardir))
from input_set import InputSet
from input_set import InputSet, MGInputSet
from openmc.statepoint import StatePoint
from openmc.executor import Executor
import openmc.particle_restart as pr
@ -243,12 +243,14 @@ class ParticleRestartTestHarness(TestHarness):
class PyAPITestHarness(TestHarness):
def __init__(self, statepoint_name, tallies_present=False):
def __init__(self, statepoint_name, tallies_present=False, mg=False):
super(PyAPITestHarness, self).__init__(statepoint_name, tallies_present)
self.parser.add_option('--build-inputs', dest='build_only',
action='store_true', default=False)
self._input_set = InputSet()
if mg:
self._input_set = MGInputSet()
else:
self._input_set = InputSet()
def main(self):
"""Accept commandline arguments and either run or update tests."""
(self._opts, self._args) = self.parser.parse_args()