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

Merge branch 'develop' into scopatz-six

Browse source
This commit is contained in:
Paul Romano 2016-10-18 22:49:56 -05:00
commit 0454e66a11
55 changed files with 4195 additions and 16187 deletions
Download patch file Download diff file Expand all files Collapse all files

408
docs/source/io_formats/mgxs_library.rst
View file

@ -6,21 +6,19 @@ 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,
``cross_sections.xml`` file contains necessary meta-data for each dataset,
including the name and a file system location where the complete library
can be found. In multi-group mode, this ``mgxs.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 ``mgxs.xml``
file.
can be found. In multi-group mode, the multi-group meta-data and the
nuclear data itself is contained within an ``mgxs.h5`` file. This portion of
the manual describes the format of the multi-group data library required
to be used in the ``mgxs.h5`` 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.
The multi-group library is provided in the HDF5_ 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/
.. _HDF5: http://www.hdfgroup.org/HDF5/
.. _mgxs_lib_spec:
@ -28,277 +26,139 @@ materials.
MGXS Library Specification
--------------------------
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:
:Attributes: - **groups** (*int*) -- Number of energy groups
- **group structure** (*double[]*) -- Monotonically increasing
list of group boundaries, in units of MeV. The length of this
array should be the number of groups plus 1.
``<groups>`` Element
--------------------
**/<library name>/**
The ``<groups>`` element has no attributes and simply provides the number of
energy groups contained within the library.
The data within <library name> contains the temperature-dependent multi-group
data for the nuclide or material that it represents.
*Default*: None, this must be provided.
:Attributes: - **atomic_weight_ratio** (*double*) -- The atomic weight ratio
(optional, i.e. it is not meaningful for material-wise data).
- **fissionable** (*bool*) -- Whether the dataset is fissionable
(True) or not (False).
- **representation** (*char[]*) -- 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 a
similar representation) and thus are angle-independent, or if the
data was generated with angular dependent fluxes and thus the
data is angle-dependent. Valid values are either "isotropic" or
"angle".
- **num_azimuthal** (*int*) -- Number of equal width angular bins
that the azimuthal angular domain is subdivided if the
`representation` attribute is "angle". This parameter is
ignored otherwise.
- **num_polar** (*int*) -- Number of equal width angular bins
that the polar angular domain is subdivided if the
`representation` attribute is "angle". This parameter is
ignored otherwise.
- **scatter_format** (*char[]*) -- The representation of the
scattering angular distribution. The options are either
"legendre", "histogram", or "tabular". If not provided, the
default of "legendre" will be assumed.
- **order** (*int*) -- Either the Legendre order, number of bins,
or number of points (depending on the value of `scatter_format`)
used to describe the angular distribution associated with each
group-to-group transfer probability.
- **scatter_shape** (*char[]*) -- The shape of the provided
scatter and multiplicity matrix. The values provided are strings
describing the ordering the scattering array is provided in
row-major (i.e., C/C++ and Python) indexing. Valid values are
"[Order][G][G']" or "[Order][G'][G]" where "G'" denotes the
secondary/outgoing energy groups, "G" denotes the incoming
energy groups, and "Order" is the angular distribution index.
This value is not required; if not the default value of
"[Order][G][G']" will be assumed.
.. _group_structure:
**/<library name>/kTs/**
``<group_structure>`` Element
-----------------------------
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin), rounded to the nearest integer
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:
**/<library name>/<TTT>K/**
.. code-block:: xml
Temperature-dependent data, provided for temperature <TTT>K.
<group_structure> 0.0 0.625E-6 20.0 </group_structure>
:Datasets: - **total** (*double[]* or *double[][][]*) -- Total cross section.
This is a 1-D vector if `representation` is "isotropic", or a 3-D
vector if `representation` is "angle" with dimensions of
[groups][azimuthal][polar].
- **absorption** (*double[]* or *double[][][]*) -- Absorption
cross section.
This is a 1-D vector if `representation` is "isotropic", or a 3-D
vector if `representation` is "angle" with dimensions of
[groups][azimuthal][polar].
- **fission** (*double[]* or *double[][][]*) -- Fission
cross section.
This is a 1-D vector if `representation` is "isotropic", or a 3-D
vector if `representation` is "angle" with dimensions of
[groups][azimuthal][polar]. This is only required if the dataset
is fissionable and fission-tallies are expected to be used.
- **kappa-fission** (*double[]* or *double[][][]*) -- Kappa-Fission
(energy-release from fission) cross section.
This is a 1-D vector if `representation` is "isotropic", or a 3-D
vector if `representation` is "angle" with dimensions of
[groups][azimuthal][polar]. This is only required if the dataset
is fissionable and fission-tallies are expected to be used.
- **chi** (*double[]* or *double[][][]*) -- Fission neutron energy
spectra.
This is a 1-D vector if `representation` is "isotropic", or a 3-D
vector if `representation` is "angle" with dimensions of
[groups][azimuthal][polar]. This is only required if the dataset
is fissionable and fission-tallies are expected to be used.
- **nu-fission** (*double[]* to *double[][][][]*) -- Nu-Fission
cross section.
If **chi** is provided, then `nu-fission` has the same
dimensionality as `fission`. If **chi** is not provided, then
the `nu-fission` data must represent the fission neutron energy
spectra as well and thus will have one additional dimension
for the outgoing energy group. In this case, `nu-fission` has the
same dimensionality as `multiplicity_matrix`.
- **inverse_velocities** (*double[]*) -- Average inverse velocity
for each of the groups in the library. This dataset is optional.
*Default*: None, this must be provided.
**/<library name>/<TTT>K/scatter_data/**
.. _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 can allow for a decrease in 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.
Data specific to neutron scattering for the temperature <TTT>K
:Datasets: - **g_min** (*int[]* or *int[][][]*) --
Minimum (most energetic) groups with non-zero values of
the scattering matrix provided. If `scatter_shape` is
"[Order][G][G']" then `g_min` will describe the minimum values
of "G'" for each "G"; if `scatter_shape` is "[Order][G'][G]"
then `g_min` will describe the minimum values of "G" for each "G'".
These group numbers use the standard
ordering where the fastest neutron energy group is group 1 while
the slowest neutron energy group is group G.
The dimensionality of `g_min` is:
`g_min[g]`, or `g_min[num_polar][num_azimuthal][g]`.
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle".
- **g_max** (*int[]* or *int[][][]*) --
Similar to `g_min`, except this dataset describes the maximum
(least energetic) groups with non-zero values of
the scattering matrix.
- **scatter_matrix** (*double[]*) -- Flattened representation of the
scattering moment matrices. The pre-flattened array corresponds to
the shape provied in `scatter_shape`, but if `representation` is
"angle" the dimensionality in `scatter_shape` is prepended by
"[num_polar][num_azimuthal]" dimensions. The right-most energy
group dimension will only include the entries between `g_min` and
`g_max`.
dimension has a dimensionality of `g_min` to `g_max`.
- **multiplicity_matrix** (*double[]*) -- Flattened representation of
the scattering moment matrices. This dataset provides the code with
a scaling factor to account for neutrons being produced in (n,xn)
reactions. This is assumed isotropic and therefore is not repeated
for every Legendre moment or histogram/tabular bin. This dataset is
optional, if it is not provided no multiplication (i.e., values of
1.0) will be assumed.
The pre-flattened array is shapes consistent with `scatter_matrix`
except the "[Order]" dimension in `scatter_shape` is ignored since
this data is assumed isotropic.

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

File diff suppressed because one or more lines are too long

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

@ -121,10 +121,13 @@ multi-group mode.
``<cutoff>`` Element
--------------------
The ``<cutoff>`` element indicates the weight cutoff used below which particles
undergo Russian roulette. Surviving particles are assigned a user-determined
weight. Note that weight cutoffs and Russian rouletting are not turned on by
default. This element has the following attributes/sub-elements:
The ``<cutoff>`` element indicates two kinds of cutoffs. The first is the weight
cutoff used below which particles undergo Russian roulette. Surviving particles
are assigned a user-determined weight. Note that weight cutoffs and Russian
rouletting are not turned on by default. The second is the energy cutoff which
is used to kill particles under certain energy. The energy cutoff should not be
used unless you know particles under the energy are of no importance to results
you care. This element has the following attributes/sub-elements:
:weight:
The weight below which particles undergo Russian roulette.
@ -137,6 +140,11 @@ default. This element has the following attributes/sub-elements:
*Default*: 1.0
:energy:
The energy under which particles will be killed.
*Default*: 0.0
.. _eigenvalue:
``<eigenvalue>`` Element
@ -709,6 +717,36 @@ survival biasing, otherwise known as implicit capture or absorption.
*Default*: false
.. _tabular_legendre:
``<tabular_legendre>`` Element
---------------------------------
The optional ``<tabular_legendre>`` element specifies how the multi-group
Legendre scattering kernel is represented if encountered in a multi-group
problem. 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 can
allow for a decrease in 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 of a
Legendre scattering expansion to the tabular format should be performed or
not. A value of “true” means the conversion should be performed, “false”
means it will 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
.. note:: This element is only used in the multi-group :ref:`energy_mode`.
.. _temperature_default:
``<temperature_default>`` Element
@ -872,6 +910,19 @@ displayed. This element takes the following attributes:
*Default*: 5
``<create_fission_neutrons>`` Element
-------------------------------------
The ``<create_fission_neutrons>`` element indicates whether fission neutrons
should be created or not. If this element is set to "true", fission neutrons
will be created; otherwise the fission is treated as capture and no fission
neutron will be created. Note that this option is only applied to fixed source
calculation. For eigenvalue calculation, fission will always be treated as real
fission.
*Default*: true
``<volume_calc>`` Element
-------------------------

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

@ -21,44 +21,47 @@ groups = openmc.mgxs.EnergyGroups(group_edges=[1E-11, 0.0635E-6, 10.0E-6,
# Instantiate the 7-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2', groups)
uo2_xsdata.order = 0
uo2_xsdata.total = [0.1779492, 0.3298048, 0.4803882, 0.5543674,
0.3118013, 0.3951678, 0.5644058]
uo2_xsdata.absorption = [8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
3.0020E-02, 1.1126E-01, 2.8278E-01]
uo2_xsdata.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.fission = [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 = [2.005998E-02, 2.027303E-03, 1.570599E-02,
4.518301E-02, 4.334208E-02, 2.020901E-01,
5.257105E-01]
uo2_xsdata.chi = [5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07,
0.0000E+00, 0.0000E+00, 0.0000E+00]
uo2_xsdata.set_total(
[0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678,
0.5644058])
uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
3.0020E-02, 1.1126E-01, 2.8278E-01])
uo2_xsdata.set_scatter_matrix(
[[[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.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03,
1.85648E-02, 1.78084E-02, 8.30348E-02,
2.16004E-01])
uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02,
4.518301E-02, 4.334208E-02, 2.020901E-01,
5.257105E-01])
uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00,
0.0000E+00, 0.0000E+00])
h2o_xsdata = openmc.XSdata('LWTR', groups)
h2o_xsdata.order = 0
h2o_xsdata.total = [0.15920605, 0.412969593, 0.59030986, 0.58435,
0.718, 1.2544497, 2.650379]
h2o_xsdata.absorption = [6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02]
h2o_xsdata.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.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
0.718, 1.2544497, 2.650379])
h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02])
h2o_xsdata.set_scatter_matrix(
[[[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]]])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
mg_cross_sections_file.export_to_xml()
mg_cross_sections_file.export_to_hdf5()
###############################################################################
@ -129,7 +132,7 @@ geometry.export_to_xml()
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings_file = openmc.Settings()
settings_file.energy_mode = "multi-group"
settings_file.cross_sections = "./mgxs.xml"
settings_file.cross_sections = "./mgxs.h5"
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles

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

@ -1,383 +0,0 @@
<?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</name>
<alias>UO2</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 -->
<kappa_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</kappa_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</name>
<alias>MOX1</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 -->
<kappa_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</kappa_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</name>
<alias>MOX2</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 -->
<kappa_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</kappa_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</name>
<alias>MOX3</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 -->
<kappa_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</kappa_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</name>
<alias>FC</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 -->
<kappa_fission>
1.0 1.0 1.0 1.0 1.0 1.0 1.0
</kappa_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</name>
<alias>GT</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</name>
<alias>LWTR</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</name>
<alias>CR</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>

BIN
examples/xml/pincell_multigroup/mgxs.h5 Normal file
View file

Binary file not shown.

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

@ -28,13 +28,12 @@
</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>
<cross_sections>./mgxs.h5</cross_sections>
</settings>

19
openmc/mesh.py
View file

@ -284,10 +284,8 @@ class Mesh(object):
openmc.XPlane(x0=self.upper_right[0],
boundary_type=bc[1])]
if len(self.dimension) == 1:
yplanes = [openmc.YPlane(y0=np.finfo(np.float).min,
boundary_type='reflective'),
openmc.YPlane(y0=np.finfo(np.float).max,
boundary_type='reflective')]
yplanes = [openmc.YPlane(y0=-1e10, boundary_type='reflective'),
openmc.YPlane(y0=1e10, boundary_type='reflective')]
else:
yplanes = [openmc.YPlane(y0=self.lower_left[1],
boundary_type=bc[2]),
@ -295,10 +293,15 @@ class Mesh(object):
boundary_type=bc[3])]
if len(self.dimension) <= 2:
zplanes = [openmc.ZPlane(z0=np.finfo(np.float).min,
boundary_type='reflective'),
openmc.ZPlane(z0=np.finfo(np.float).max,
boundary_type='reflective')]
# Would prefer to have the z ranges be the max supported float, but
# these values are apparently different between python and Fortran.
# Choosing a safe and sane default.
# Values of +/-1e10 are used here as there seems to be an
# inconsistency between what numpy uses as the max float and what
# Fortran expects for a real(8), so this avoids code complication
# and achieves the same goal.
zplanes = [openmc.ZPlane(z0=-1e10, boundary_type='reflective'),
openmc.ZPlane(z0=1e10, boundary_type='reflective')]
else:
zplanes = [openmc.ZPlane(z0=self.lower_left[2],
boundary_type=bc[4]),

136
openmc/mgxs/library.py
View file

@ -819,11 +819,13 @@ class Library(object):
return pickle.load(open(full_filename, 'rb'))
def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
order=None, tabular_legendre=None, tabular_points=33,
subdomain=None):
order=None, subdomain=None):
"""Generates an openmc.XSdata object describing a multi-group cross section
data set for eventual combination in to an openmc.MGXSLibrary object
(i.e., the library).
dataset for writing to an openmc.MGXSLibrary object.
Note that this method does not build an XSdata
object with nested temperature tables. The temperature of each
XSdata object will be left at the default value of 300K.
Parameters
----------
@ -842,18 +844,6 @@ class Library(object):
Scattering order for this data entry. Default is None,
which will set the XSdata object to use the order of the
Library.
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
subdomain : iterable of int
This parameter is not used unless using a mesh domain. In that
case, the subdomain is an [i,j,k] index (1-based indexing) of the
@ -864,7 +854,7 @@ class Library(object):
Returns
-------
xsdata : openmc.XSdata
Multi-Group Cross Section data set object.
Multi-Group Cross Section dataset object.
Raises
------
@ -887,10 +877,6 @@ class Library(object):
if order is not None:
cv.check_greater_than('order', order, 0, equality=True)
cv.check_less_than('order', order, 10, equality=True)
cv.check_type('tabular_legendre', tabular_legendre,
(type(None), bool))
if tabular_points is not None:
cv.check_greater_than('tabular_points', tabular_points, 1)
if subdomain is not None:
cv.check_iterable_type('subdomain', subdomain, Integral,
max_depth=3)
@ -907,7 +893,7 @@ class Library(object):
# Build & add metadata to XSdata object
name = xsdata_name
if nuclide is not 'total':
if nuclide != 'total':
name += '_' + nuclide
xsdata = openmc.XSdata(name, self.energy_groups)
@ -919,14 +905,12 @@ class Library(object):
# the provided order or the Library's order.
xsdata.order = min(order, self.legendre_order)
# Set the tabular_legendre option if needed
if tabular_legendre is not None:
xsdata.tabular_legendre = {'enable': tabular_legendre,
'num_points': tabular_points}
# Right now only 'legendre' data and isotropic weighting is supported
self.scatter_format = 'legendre'
self.representation = 'isotropic'
if nuclide is not 'total':
xsdata.zaid = self._nuclides[nuclide][0]
xsdata.awr = self._nuclides[nuclide][1]
if nuclide != 'total':
xsdata.atomic_weight_ratio = self._nuclides[nuclide][1]
if subdomain is None:
subdomain = 'all'
@ -937,7 +921,7 @@ class Library(object):
if 'nu-transport' in self.mgxs_types and self.correction == 'P0':
mymgxs = self.get_mgxs(domain, 'nu-transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomains=subdomain)
subdomain=subdomain)
elif 'total' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'total')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
@ -976,49 +960,55 @@ class Library(object):
# If multiplicity matrix is available, prefer that
if 'multiplicity matrix' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'multiplicity matrix')
xsdata.set_multiplicity_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
xsdata.set_multiplicity_matrix_mgxs(mymgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
# multiplicity will fall back to using scatter and nu-scatter
elif ((('scatter matrix' in self.mgxs_types) and
('nu-scatter matrix' in self.mgxs_types))):
scatt_mgxs = self.get_mgxs(domain, 'scatter matrix')
nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
xsdata.set_multiplicity_mgxs(nuscatt_mgxs, scatt_mgxs,
xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
xsdata.set_multiplicity_matrix_mgxs(nuscatt_mgxs, scatt_mgxs,
xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
using_multiplicity = True
else:
using_multiplicity = False
if using_multiplicity:
nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
xsdata.set_scatter_mgxs(nuscatt_mgxs, xs_type=xs_type,
nuclide=[nuclide], subdomain=subdomain)
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
else:
if 'nu-scatter matrix' in self.mgxs_types:
nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
xsdata.set_scatter_mgxs(nuscatt_mgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
xsdata.set_scatter_matrix_mgxs(nuscatt_mgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
# Since we are not using multiplicity, then
# scattering multiplication (nu-scatter) must be
# accounted for approximately by using an adjusted
# absorption cross section.
if 'total' in self.mgxs_types:
xsdata.absorption = \
np.subtract(xsdata.total,
np.sum(xsdata.scatter[0, :, :], axis=1))
if 'total' in self.mgxs_types or 'transport' in self.mgxs_types:
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
return xsdata
def create_mg_library(self, xs_type='macro', xsdata_names=None,
tabular_legendre=None, tabular_points=33):
def create_mg_library(self, xs_type='macro', xsdata_names=None):
"""Creates an openmc.MGXSLibrary object to contain the MGXS data for the
Multi-Group mode of OpenMC.
Note that this library will not make use of nested temperature tables.
Every dataset in the library will be treated as if it was at the same
default temperature.
Parameters
----------
xs_type: {'macro', 'micro'}
@ -1028,18 +1018,6 @@ class Library(object):
xsdata_names : Iterable of str
List of names to apply to the "xsdata" entries in the
resultant mgxs data file. Defaults to 'set1', 'set2', ...
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
Returns
-------
@ -1091,8 +1069,6 @@ class Library(object):
# Create XSdata and Macroscopic for this domain
xsdata = self.get_xsdata(domain, xsdata_name,
tabular_legendre=tabular_legendre,
tabular_points=tabular_points,
subdomain=subdomain)
mgxs_file.add_xsdata(xsdata)
i += 1
@ -1110,45 +1086,34 @@ class Library(object):
xsdata_name = 'set' + str(i + 1)
else:
xsdata_name = xsdata_names[i]
if nuclide is not 'total':
if nuclide != 'total':
xsdata_name += '_' + nuclide
xsdata = self.get_xsdata(domain, xsdata_name,
nuclide=nuclide, xs_type=xs_type,
tabular_legendre=tabular_legendre,
tabular_points=tabular_points)
nuclide=nuclide, xs_type=xs_type)
mgxs_file.add_xsdata(xsdata)
return mgxs_file
def create_mg_mode(self, xsdata_names=None, tabular_legendre=None,
tabular_points=33, bc=['reflective'] * 6):
def create_mg_mode(self, xsdata_names=None, bc=['reflective'] * 6):
"""Creates an openmc.MGXSLibrary object to contain the MGXS data for the
Multi-Group mode of OpenMC as well as the associated openmc.Materials
and openmc.Geometry objects. The created Geometry is the same as that
used to generate the MGXS data, with the only differences being
modifications to point to newly-created Materials which point to the
multi-group data. This method only creates a macroscopic
MGXS Library even if nuclidic tallies are specified in the Library.
and openmc.Geometry objects.
The created Geometry is the same as that used to generate the MGXS
data, with the only differences being modifications to point to
newly-created Materials which point to the multi-group data. This
method only creates a macroscopic MGXS Library even if nuclidic tallies
are specified in the Library. Note that this library will not make
use of nested temperature tables. Every dataset in the library will
be treated as if it was at the same default temperature.
Parameters
----------
xsdata_names : Iterable of str
List of names to apply to the "xsdata" entries in the
resultant mgxs data file. Defaults to 'set1', 'set2', ...
tabular_legendre : None or bool
Flag to denote whether or not the Legendre expansion of the
scattering angular distribution is to be converted to a tabular
representation by OpenMC. A value of `True` means that it is to be
converted while a value of `False` means that it will not be.
Defaults to `None` which leaves the default behavior of OpenMC in
place (the distribution is converted to a tabular representation).
tabular_points : int
This parameter is not used unless the ``tabular_legendre``
parameter is set to `True`. In this case, this parameter sets the
number of equally-spaced points in the domain of [-1,1] to be used
in building the tabular distribution. Default is `33`.
bc : iterable of {'reflective', 'periodic', 'transmission', or 'vacuum'}
Boundary conditions for each of the four faces of a rectangle
(if applying to a 2D mesh) or six faces of a parallelepiped
@ -1195,8 +1160,7 @@ class Library(object):
cv.check_length("domains", self.domains, 1, 1)
# Get the MGXS File Data
mgxs_file = self.create_mg_library('macro', xsdata_names,
tabular_legendre, tabular_points)
mgxs_file = self.create_mg_library('macro', xsdata_names)
# Now move on the creating the geometry and assigning materials
if self.domain_type == 'mesh':

1150
openmc/mgxs_library.py
View file

File diff suppressed because it is too large Load diff

126
openmc/settings.py
View file

@ -87,10 +87,14 @@ class Settings(object):
Seed for the linear congruential pseudorandom number generator
survival_biasing : bool
Indicate whether survival biasing is to be used
weight : float
Weight cutoff below which particle undergo Russian roulette
weight_avg : float
Weight assigned to particles that are not killed after Russian roulette
cutoff : dict
Dictionary defining weight cutoff and energy cutoff. The dictionary may
have three keys, 'weight', 'weight_avg' and 'energy'. Value for 'weight'
should be a float indicating weight cutoff below which particle undergo
Russian roulette. Value for 'weight_avg' should be a float indicating
weight assigned to particles that are not killed after Russian
roulette. Value of energy should be a float indicating energy in MeV
below which particle will be killed.
entropy_dimension : tuple or list
Number of Shannon entropy mesh cells in the x, y, and z directions,
respectively
@ -98,6 +102,13 @@ class Settings(object):
Coordinates of the lower-left point of the Shannon entropy mesh
entropy_upper_right : tuple or list
Coordinates of the upper-right point of the Shannon entropy mesh
tabular_legendre : dict
Determines if a multi-group scattering moment kernel expanded via
Legendre polynomials is to be converted to a tabular distribution or
not. Accepted keys are 'enable' and 'num_points'. The value for
'enable' is a bool stating whether the conversion to tabular is
performed; the value for 'num_points' sets the number of points to use
in the tabular distribution, should 'enable' be True.
temperature : dict
Defines a default temperature and method for treating intermediate
temperatures at which nuclear data doesn't exist. Accepted keys are
@ -139,6 +150,8 @@ class Settings(object):
The elastic scattering model to use for resonant isotopes
volume_calculations : VolumeCalculation or iterable of VolumeCalculation
Stochastic volume calculation specifications
create_fission_neutrons : bool
Indicate whether fission neutrons should be created or not.
"""
@ -198,11 +211,12 @@ class Settings(object):
self._trace = None
self._track = None
self._tabular_legendre = {}
self._temperature = {}
# Cutoff subelement
self._weight = None
self._weight_avg = None
self._cutoff = None
# Uniform fission source subelement
self._ufs_dimension = 1
@ -225,6 +239,8 @@ class Settings(object):
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'volume calculations')
self._create_fission_neutrons = None
@property
def run_mode(self):
return self._run_mode
@ -361,6 +377,10 @@ class Settings(object):
def verbosity(self):
return self._verbosity
@property
def tabular_legendre(self):
return self._tabular_legendre
@property
def temperature(self):
return self._temperature
@ -374,12 +394,8 @@ class Settings(object):
return self._track
@property
def weight(self):
return self._weight
@property
def weight_avg(self):
return self._weight_avg
def cutoff(self):
return self._cutoff
@property
def ufs_dimension(self):
@ -425,6 +441,10 @@ class Settings(object):
def volume_calculations(self):
return self._volume_calculations
@property
def create_fission_neutrons(self):
return self._create_fission_neutrons
@run_mode.setter
def run_mode(self, run_mode):
if run_mode not in ['eigenvalue', 'fixed source']:
@ -612,17 +632,29 @@ class Settings(object):
cv.check_type('survival biasing', survival_biasing, bool)
self._survival_biasing = survival_biasing
@weight.setter
def weight(self, weight):
cv.check_type('weight cutoff', weight, Real)
cv.check_greater_than('weight cutoff', weight, 0.0)
self._weight = weight
@cutoff.setter
def cutoff(self, cutoff):
if not isinstance(cutoff, Mapping):
msg = 'Unable to set cutoff from "{0}" which is not a '\
' Python dictionary'.format(cutoff)
raise ValueError(msg)
for key in cutoff:
if key == 'weight':
cv.check_type('weight cutoff', cutoff['weight'], Real)
cv.check_greater_than('weight cutoff', cutoff['weight'], 0.0)
elif key == 'weight_avg':
cv.check_type('average survival weight', cutoff['weight_avg'],
Real)
cv.check_greater_than('average survival weight',
cutoff['weight_avg'], 0.0)
elif key == 'energy':
cv.check_type('energy cutoff', cutoff['energy'], Real)
cv.check_greater_than('energy cutoff', cutoff['energy'], 0.0)
else:
msg = 'Unable to set cutoff to "{0}" which is unsupported by '\
'OpenMC'.format(key)
@weight_avg.setter
def weight_avg(self, weight_avg):
cv.check_type('average survival weight', weight_avg, Real)
cv.check_greater_than('average survival weight', weight_avg, 0.0)
self._weight_avg = weight_avg
self._cutoff = cutoff
@entropy_dimension.setter
def entropy_dimension(self, dimension):
@ -666,6 +698,19 @@ class Settings(object):
cv.check_type('no reduction option', no_reduce, bool)
self._no_reduce = no_reduce
@tabular_legendre.setter
def tabular_legendre(self, tabular_legendre):
cv.check_type('tabular_legendre settings', tabular_legendre, Mapping)
for key, value in tabular_legendre.items():
cv.check_value('tabular_legendre key', key,
['enable', 'num_points'])
if key == 'enable':
cv.check_type('enable tabular_legendre', value, bool)
elif key == 'num_points':
cv.check_type('num_points tabular_legendre', value, Integral)
cv.check_greater_than('num_points tabular_legendre', value, 0)
self._tabular_legendre = tabular_legendre
@temperature.setter
def temperature(self, temperature):
cv.check_type('temperature settings', temperature, Mapping)
@ -824,6 +869,12 @@ class Settings(object):
self._volume_calculations = cv.CheckedList(
VolumeCalculation, 'stochastic volume calculations', vol_calcs)
@create_fission_neutrons.setter
def create_fission_neutrons(self, create_fission_neutrons):
cv.check_type('Whether create fission neutrons',
create_fission_neutrons, bool)
self._create_fission_neutrons = create_fission_neutrons
def _create_run_mode_subelement(self):
if self.run_mode == 'eigenvalue':
@ -988,14 +1039,19 @@ class Settings(object):
element.text = str(self._survival_biasing).lower()
def _create_cutoff_subelement(self):
if self._weight is not None:
if self._cutoff is not None:
element = ET.SubElement(self._settings_file, "cutoff")
if 'weight' in self._cutoff:
subelement = ET.SubElement(element, "weight")
subelement.text = str(self._cutoff['weight'])
subelement = ET.SubElement(element, "weight")
subelement.text = str(self._weight)
if 'weight_avg' in self._cutoff:
subelement = ET.SubElement(element, "weight_avg")
subelement.text = str(self._cutoff['weight_avg'])
subelement = ET.SubElement(element, "weight_avg")
subelement.text = str(self._weight_avg)
if 'energy' in self._cutoff:
subelement = ET.SubElement(element, "energy")
subelement.text = str(self._cutoff['energy'])
def _create_entropy_subelement(self):
if self._entropy_lower_left is not None and \
@ -1049,6 +1105,15 @@ class Settings(object):
element = ET.SubElement(self._settings_file, "no_reduce")
element.text = str(self._no_reduce).lower()
def _create_tabular_legendre_subelements(self):
if self.tabular_legendre:
element = ET.SubElement(self._settings_file, "tabular_legendre")
subelement = ET.SubElement(element, "enable")
subelement.text = str(self._tabular_legendre['enable']).lower()
if 'num_points' in self._tabular_legendre:
subelement = ET.SubElement(element, "num_points")
subelement.text = str(self._tabular_legendre['num_points'])
def _create_temperature_subelements(self):
if self.temperature:
for key, value in sorted(self.temperature.items()):
@ -1119,6 +1184,11 @@ class Settings(object):
for r in self.resonance_scattering:
elem.append(r.to_xml_element())
def _create_create_fission_neutrons_subelement(self):
if self._create_fission_neutrons is not None:
elem = ET.SubElement(self._settings_file, "create_fission_neutrons")
elem.text = str(self._create_fission_neutrons).lower()
def export_to_xml(self, path='settings.xml'):
"""Export simulation settings to an XML file.
@ -1155,6 +1225,7 @@ class Settings(object):
self._create_no_reduce_subelement()
self._create_threads_subelement()
self._create_verbosity_subelement()
self._create_tabular_legendre_subelements()
self._create_temperature_subelements()
self._create_trace_subelement()
self._create_track_subelement()
@ -1162,6 +1233,7 @@ class Settings(object):
self._create_dd_subelement()
self._create_resonance_scattering_subelement()
self._create_volume_calcs_subelement()
self._create_create_fission_neutrons_subelement()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._settings_file)

2
scripts/openmc-ace-to-hdf5
View file

@ -108,7 +108,7 @@ elif args.xsdata is not None:
for line in xsdata:
words = line.split()
if len(words) >= 9:
path = os.path.join(os.path.dirname(args.xsdata, words[8]))
path = os.path.join(os.path.dirname(args.xsdata), words[8])
if path not in ace_libraries:
ace_libraries.append(path)

229
scripts/openmc-update-mgxs Executable file
View file

@ -0,0 +1,229 @@
#!/usr/bin/env python
"""Update OpenMC's deprecated multi-group cross section XML files to the latest
HDF5-based format.
Usage information can be obtained by running 'openmc-update-mgxs --help':
usage: openmc-update-mgxs [-h] in out
Update mgxs.xml files to the latest format. This will remove 'outside'
attributes/elements from lattices and replace them with 'outer' attributes. For
'cell' elements, any 'surfaces' attributes/elements will be renamed
'region'. Note that this script will not delete the given files; it will append
'.original' to the given files and write new ones.
positional arguments:
in Input mgxs xml file
out Output mgxs hdf5 file
optional arguments:
-h, --help show this help message and exit
"""
from __future__ import print_function
import os
from shutil import move
import warnings
import xml.etree.ElementTree as ET
import argparse
import h5py
import numpy as np
import openmc.mgxs_library
description = """\
Update OpenMC's deprecated multi-group cross section XML files to the latest
HDF5-based format."""
def parse_args():
"""Read the input files from the commandline."""
# Create argument parser
parser = argparse.ArgumentParser(description=description,
formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument('-i', '--input', type=argparse.FileType('r'),
help='input XML file')
parser.add_argument('-o', '--output', nargs='?', default='',
help='output file, in HDF5 format')
args = vars(parser.parse_args())
if args['output'] == '':
filename = args['input'].name
extension = filenameos.path.splitext()
if extension == '.xml':
filename = filename[:filename.rfind('.')] + '.h5'
args['output'] = filename
# Parse and return commandline arguments.
return args
def get_data(element, entry):
value = element.find(entry)
if value is not None:
value = value.text.strip()
else:
if entry in element.attrib:
value = element.attrib[entry].strip()
else:
value = None
return value
if __name__ == '__main__':
args = parse_args()
# Parse the XML data.
tree = ET.parse(args['input'])
root = tree.getroot()
# Get old metadata
temp = tree.find('group_structure').text.strip()
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse_velocities')
if temp is not None:
temp = temp.text.strip()
temp = np.array(temp.split())
inverse_velocities = temp.astype(np.float)
else:
inverse_velocities = None
xsd = []
names = []
# Now move on to the cross section data itself
for xsdata_elem in root.iter('xsdata'):
name = get_data(xsdata_elem, 'name')
temperature = get_data(xsdata_elem, 'kT')
if temperature is not None:
temperature = \
float(temperature) / openmc.data.K_BOLTZMANN
else:
temperature = 294.
temperatures = [temperature]
awr = get_data(xsdata_elem, 'awr')
if awr is not None:
awr = float(awr)
representation = get_data(xsdata_elem, 'representation')
if representation is None:
representation = 'isotropic'
if representation == 'angle':
n_azi = int(get_data(xsdata_elem, 'num_azimuthal'))
n_pol = int(get_data(xsdata_elem, 'num_polar'))
scatter_format = get_data(xsdata_elem, 'scatt_type')
if scatter_format is None:
scatter_format = 'legendre'
order = int(get_data(xsdata_elem, 'order'))
tab_leg = get_data(xsdata_elem, 'tabular_legendre')
if tab_leg is not None:
warnings.Warning('The tabular_legendre option has moved to the '
'settings.xml file and must be added manually')
# Either add the data to a previously existing xsdata (if it is
# for the same 'name' but a different temperature), or create a
# new one.
try:
# It is in our list, so store that entry
i = names.index(name)
except ValueError:
# It is not in our list, so add it
i = -1
xsd.append(openmc.XSdata(name, energy_groups,
temperatures=temperatures,
representation=representation))
if awr is not None:
xsd[-1].atomic_weight_ratio = awr
if representation == 'angle':
xsd[-1].num_azimuthal = n_azi
xsd[-1].num_polar = n_pol
xsd[-1].scatter_format = scatter_format
xsd[-1].order = order
names.append(name)
if scatter_format == 'legendre':
order_dim = order + 1
else:
order_dim = order
if i != -1:
xsd[i].add_temperature(temperature)
temp = get_data(xsdata_elem, 'total')
if temp is not None:
temp = np.array(temp.split(), dtype=float)
total = temp.astype(np.float)
total.shape = xsd[i].vector_shape
xsd[i].set_total(total, temperature)
if inverse_velocities is not None:
xsd[i].set_inverse_velocities(inverse_velocities, temperature)
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())
absorption = temp.astype(np.float)
absorption.shape = xsd[i].vector_shape
xsd[i].set_absorption(absorption, temperature)
temp = get_data(xsdata_elem, 'scatter')
temp = np.array(temp.split())
scatter = temp.astype(np.float)
scatter.shape = xsd[i].pn_matrix_shape
xsd[i].set_scatter_matrix(scatter, temperature)
temp = get_data(xsdata_elem, 'multiplicity')
if temp is not None:
temp = np.array(temp.split())
multiplicity = temp.astype(np.float)
multiplicity.shape = xsd[i].matrix_shape
xsd[i].set_multiplicity_matrix(multiplicity, temperature)
temp = get_data(xsdata_elem, 'fission')
if temp is not None:
temp = np.array(temp.split())
fission = temp.astype(np.float)
fission.shape = xsd[i].vector_shape
xsd[i].set_fission(fission, temperature)
temp = get_data(xsdata_elem, 'kappa_fission')
if temp is not None:
temp = np.array(temp.split())
kappa_fission = temp.astype(np.float)
kappa_fission.shape = xsd[i].vector_shape
xsd[i].set_kappa_fission(kappa_fission, temperature)
temp = get_data(xsdata_elem, 'chi')
if temp is not None:
temp = np.array(temp.split())
chi = temp.astype(np.float)
chi.shape = xsd[i].vector_shape
xsd[i].set_chi(chi, temperature)
else:
chi = None
temp = get_data(xsdata_elem, 'nu_fission')
if temp is not None:
temp = np.array(temp.split())
nu_fission = temp.astype(np.float)
if chi is not None:
nu_fission.shape = xsd[i].vector_shape
else:
nu_fission.shape = xsd[i].matrix_shape
xsd[i].set_nu_fission(nu_fission, temperature)
# Build library as we go, but first we have enough to initialize it
lib = openmc.MGXSLibrary(energy_groups)
lib.add_xsdatas(xsd)
lib.export_to_hdf5(args['output'])

4
setup.py Normal file → Executable file
View file

@ -1,6 +1,7 @@
#!/usr/bin/env python
import glob
import numpy as np
try:
from setuptools import setup
have_setuptools = True
@ -57,7 +58,8 @@ if have_setuptools:
# If Cython is present, add resonance reconstruction capability
if have_cython:
kwargs.update({
'ext_modules': cythonize('openmc/data/reconstruct.pyx')
'ext_modules': cythonize('openmc/data/reconstruct.pyx'),
'include_dirs': [np.get_include()]
})
setup(**kwargs)

78
src/geometry_header.F90
View file

@ -1,6 +1,12 @@
module geometry_header
use constants, only: HALF, TWO, THREE, INFINITY
use algorithm, only: find
use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, &
MATERIAL_VOID, NONE
use dict_header, only: DictCharInt, DictIntInt
use material_header, only: Material
use stl_vector, only: VectorReal
use string, only: to_lower
implicit none
@ -319,4 +325,74 @@ contains
end if
end function get_local_hex
!===============================================================================
! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table
! appears at in the model. Later, this list is used to determine the actual
! temperatures to read (which may be different if interpolation is used)
!===============================================================================
subroutine get_temperatures(cells, materials, material_dict, nuclide_dict, &
n_nucs, nuc_temps, sab_dict, n_sabs, sab_temps)
type(Cell), allocatable, intent(in) :: cells(:)
type(Material), allocatable, intent(in) :: materials(:)
type(DictIntInt), intent(in) :: material_dict
type(DictCharInt), intent(in) :: nuclide_dict
integer, intent(in) :: n_nucs
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(DictCharInt), optional, intent(in) :: sab_dict
integer, optional, intent(in) :: n_sabs
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_sab ! index in S(a,b) array
integer :: i_material
real(8) :: temperature ! temperature in Kelvin
allocate(nuc_temps(n_nucs))
if (present(n_sabs) .and. present(sab_temps)) allocate(sab_temps(n_sabs))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(cells(i) % material(j))
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in nuc_temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
! Add temperature if it hasn't already been added
if (find(nuc_temps(i_nuclide), temperature) == -1) then
call nuc_temps(i_nuclide) % push_back(temperature)
end if
end do NUC_NAMES_LOOP
if (present(sab_temps) .and. present(sab_dict) .and. &
mat % n_sab > 0) then
SAB_NAMES_LOOP: do k = 1, size(mat % sab_names)
! Get index in nuc_temps array
i_sab = sab_dict % get_key(to_lower(mat % sab_names(k)))
! Add temperature if it hasn't already been added
if (find(sab_temps(i_sab), temperature) == -1) then
call sab_temps(i_sab) % push_back(temperature)
end if
end do SAB_NAMES_LOOP
end if
end associate
end do
end do
end subroutine get_temperatures
end module geometry_header

13
src/global.F90
View file

@ -123,12 +123,15 @@ module global
! 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
! Whether or not to convert Legendres to tabulars
logical :: legendre_to_tabular = .True.
! Number of points to use in the Legendre to tabular conversion
integer :: legendre_to_tabular_points = 33
! ============================================================================
! TALLY-RELATED VARIABLES
@ -304,6 +307,7 @@ module global
logical :: survival_biasing = .false.
real(8) :: weight_cutoff = 0.25_8
real(8) :: energy_cutoff = ZERO
real(8) :: weight_survive = ONE
! ============================================================================
@ -352,6 +356,9 @@ module global
! Write out initial source
logical :: write_initial_source = .false.
! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE
logical :: create_fission_neutrons = .true.
! ============================================================================
! CMFD VARIABLES

95
src/hdf5_interface.F90
View file

@ -74,6 +74,7 @@ module hdf5_interface
module procedure read_attribute_integer_2D
module procedure read_attribute_string
module procedure read_attribute_string_1D
module procedure read_attribute_logical
end interface read_attribute
interface write_attribute
@ -84,12 +85,14 @@ module hdf5_interface
public :: write_dataset
public :: read_dataset
public :: attribute_exists
public :: write_attribute
public :: read_attribute
public :: file_create
public :: file_open
public :: file_close
public :: create_group
public :: object_exists
public :: open_group
public :: close_group
public :: open_dataset
@ -98,6 +101,7 @@ module hdf5_interface
public :: write_attribute_string
public :: get_groups
public :: get_datasets
public :: get_name
contains
@ -213,12 +217,11 @@ contains
subroutine get_groups(object_id, names)
integer(HID_T), intent(in) :: object_id
character(len=255), allocatable, intent(out) :: names(:)
character(len=150), allocatable, intent(out) :: names(:)
integer :: n_members, i, group_count, type
integer :: hdf5_err
character(len=255) :: name
character(len=150) :: name
! Get number of members in this location
call h5gn_members_f(object_id, './', n_members, hdf5_err)
@ -245,17 +248,49 @@ contains
end subroutine get_groups
!===============================================================================
! CHECK_ATTRIBUTE Checks to see if an attribute exists in the object
!===============================================================================
function attribute_exists(object_id, name) result(exists)
integer(HID_T), intent(in) :: object_id
character(*), intent(in) :: name ! name of group
logical :: exists
integer :: hdf5_err ! HDF5 error code
! Check if attribute exists
call h5aexists_by_name_f(object_id, '.', trim(name), exists, hdf5_err)
end function attribute_exists
!===============================================================================
! CHECK_GROUP Checks to see if a group exists in the object
!===============================================================================
function object_exists(object_id, name) result(exists)
integer(HID_T), intent(in) :: object_id
character(*), intent(in) :: name ! name of group
logical :: exists
integer :: hdf5_err ! HDF5 error code
! Check if group exists
call h5ltpath_valid_f(object_id, trim(name), .true., exists, hdf5_err)
end function object_exists
!===============================================================================
! GET_DATASETS Gets a list of all the datasets in a given location.
!===============================================================================
subroutine get_datasets(object_id, names)
integer(HID_T), intent(in) :: object_id
character(len=255), allocatable, intent(out) :: names(:)
character(len=150), allocatable, intent(out) :: names(:)
integer :: n_members, i, dset_count, type
integer :: hdf5_err
character(len=255) :: name
character(len=150) :: name
! Get number of members in this location
@ -283,6 +318,27 @@ contains
end subroutine get_datasets
!===============================================================================
! GET_NAME Obtains the name of the current group in group_id
!===============================================================================
function get_name(group_id, name_len_) result(name)
integer(HID_T), intent(in) :: group_id
integer(SIZE_T), optional, intent(in) :: name_len_
character(len=150) :: name ! name of group
integer(SIZE_T) :: name_len, name_file_len
integer :: hdf5_err ! HDF5 error code
if (present(name_len_)) then
name_len = name_len_
else
name_len = 150
end if
call h5iget_name_f(group_id, name, name_len, name_file_len, hdf5_err)
end function get_name
!===============================================================================
! OPEN_GROUP opens an existing HDF5 group
!===============================================================================
@ -296,7 +352,7 @@ contains
integer :: hdf5_err ! HDF5 error code
! Check if group exists
call h5ltpath_valid_f(group_id, trim(name), .true., exists, hdf5_err)
exists = object_exists(group_id, name)
! open group if it exists
if (exists) then
@ -319,7 +375,7 @@ contains
logical :: exists ! does the group exist
! Check if group exists
call h5ltpath_valid_f(group_id, trim(name), .true., exists, hdf5_err)
exists = object_exists(group_id, name)
! create group
if (exists) then
@ -357,7 +413,7 @@ contains
integer :: hdf5_err ! HDF5 error code
! Check if group exists
call h5ltpath_valid_f(group_id, trim(name), .true., exists, hdf5_err)
exists = object_exists(group_id, name)
! open group if it exists
if (exists) then
@ -2486,6 +2542,29 @@ contains
call h5tclose_f(memtype, hdf5_err)
end subroutine read_attribute_string_1D_explicit
subroutine read_attribute_logical(buffer, obj_id, name)
logical, intent(inout), target :: buffer
integer(HID_T), intent(in) :: obj_id
character(*), intent(in) :: name
integer, target :: int_buffer
integer :: hdf5_err
integer(HID_T) :: attr_id
type(c_ptr) :: f_ptr
call h5aopen_f(obj_id, trim(name), attr_id, hdf5_err)
f_ptr = c_loc(int_buffer)
call h5aread_f(attr_id, H5T_NATIVE_INTEGER, f_ptr, hdf5_err)
call h5aclose_f(attr_id, hdf5_err)
! Convert to Fortran logical
if (int_buffer == 0) then
buffer = .false.
else
buffer = .true.
end if
end subroutine read_attribute_logical
subroutine get_shape(obj_id, dims)
integer(HID_T), intent(in) :: obj_id
integer(HSIZE_T), intent(out) :: dims(:)

6
src/initialize.F90
View file

@ -113,12 +113,6 @@ contains
if (run_CE) then
! Construct log energy grid for cross-sections
call logarithmic_grid()
else
! Create material macroscopic data for MGXS
call time_read_xs%start()
call read_mgxs()
call create_macro_xs()
call time_read_xs%stop()
end if
! Allocate and setup tally stride, matching_bins, and tally maps

226
src/input_xml.F90
View file

@ -1,21 +1,20 @@
module input_xml
use hdf5
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures
use global
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
use mesh_header, only: RegularMesh
use mgxs_data, only: create_macro_xs, read_mgxs
use multipole, only: multipole_read
use output, only: write_message
use plot_header
@ -52,6 +51,17 @@ contains
call read_tallies_xml()
if (cmfd_run) call configure_cmfd()
if (.not. run_CE) then
! Create material macroscopic data for MGXS
call time_read_xs % start()
call read_mgxs()
call create_macro_xs()
call time_read_xs % stop()
end if
! Normalize atom/weight percents
if (run_mode /= MODE_PLOTTING) call normalize_ao()
end subroutine read_input_xml
!===============================================================================
@ -91,6 +101,7 @@ contains
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_keff_trigger => null()
type(Node), pointer :: node_vol => null()
type(Node), pointer :: node_tab_leg => null()
type(NodeList), pointer :: node_scat_list => null()
type(NodeList), pointer :: node_source_list => null()
type(NodeList), pointer :: node_vol_list => null()
@ -659,8 +670,15 @@ contains
! Cutoffs
if (check_for_node(doc, "cutoff")) then
call get_node_ptr(doc, "cutoff", node_cutoff)
call get_node_value(node_cutoff, "weight", weight_cutoff)
call get_node_value(node_cutoff, "weight_avg", weight_survive)
if (check_for_node(node_cutoff, "weight")) then
call get_node_value(node_cutoff, "weight", weight_cutoff)
end if
if (check_for_node(node_cutoff, "weight_avg")) then
call get_node_value(node_cutoff, "weight_avg", weight_survive)
end if
if (check_for_node(node_cutoff, "energy")) then
call get_node_value(node_cutoff, "energy", energy_cutoff)
end if
end if
! Particle trace
@ -1102,6 +1120,44 @@ contains
end select
end if
! Check whether create fission sites
if (run_mode == MODE_FIXEDSOURCE) then
if (check_for_node(doc, "create_fission_neutrons")) then
call get_node_value(doc, "create_fission_neutrons", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then
create_fission_neutrons = .true.
else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then
create_fission_neutrons = .false.
end if
end if
end if
! Check for tabular_legendre options
if (check_for_node(doc, "tabular_legendre")) then
! Get pointer to tabular_legendre node
call get_node_ptr(doc, "tabular_legendre", node_tab_leg)
! Check for enable option
if (check_for_node(node_tab_leg, "enable")) then
call get_node_value(node_tab_leg, "enable", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'false' .or. &
trim(temp_str) == '0') legendre_to_tabular = .false.
end if
! Check for the number of points
if (check_for_node(node_tab_leg, "num_points")) then
call get_node_value(node_tab_leg, "num_points", &
legendre_to_tabular_points)
if (legendre_to_tabular_points <= 1 .and. (.not. run_CE)) then
call fatal_error("The 'num_points' subelement/attribute of the &
&'tabular_legendre' element must contain a value greater than 1")
end if
end if
end if
! Close settings XML file
call close_xmldoc(doc)
@ -1366,11 +1422,7 @@ contains
! Read cell temperatures. If the temperature is not specified, set it to
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
! the temperature from the material data.
if (.not. run_CE) then
! Cell temperatures are not used for MG mode.
allocate(c % sqrtkT(1))
c % sqrtkT(1) = ZERO
else if (check_for_node(node_cell, "temperature")) then
if (check_for_node(node_cell, "temperature")) then
n = get_arraysize_double(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
@ -2048,7 +2100,7 @@ contains
if (run_CE) then
call read_ce_cross_sections_xml(libraries)
else
call read_mg_cross_sections_xml(libraries)
call read_mg_cross_sections_header(libraries)
end if
! Creating dictionary that maps the name of the material to the entry
@ -2076,18 +2128,17 @@ contains
call assign_temperatures(material_temps)
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
call get_temperatures(cells, materials, material_dict, nuclide_dict, &
n_nuclides_total, nuc_temps, sab_dict, &
n_sab_tables, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs%start()
call time_read_xs % start()
call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
call time_read_xs%stop()
call time_read_xs % stop()
end if
! Normalize atom/weight percents
if (run_mode /= MODE_PLOTTING) call normalize_ao()
! Clear dictionary
call library_dict % clear()
end subroutine read_materials
@ -3250,7 +3301,7 @@ contains
! Search through nuclides
pair_list => nuclide_dict % keys()
do while (associated(pair_list))
if (starts_with(pair_list % key, word)) then
if (trim(pair_list % key) == trim(word)) then
word = pair_list % key(1:150)
exit
end if
@ -4587,44 +4638,43 @@ contains
end subroutine read_ce_cross_sections_xml
subroutine read_mg_cross_sections_xml(libraries)
subroutine read_mg_cross_sections_header(libraries)
type(Library), allocatable, intent(out) :: libraries(:)
integer :: i ! loop index
integer :: n_libraries
logical :: file_exists ! does cross_sections.xml exist?
type(Node), pointer :: doc => null()
type(Node), pointer :: node_xsdata => null()
type(NodeList), pointer :: node_xsdata_list => null()
logical :: file_exists ! does mgxs.h5 exist?
integer(HID_T) :: file_id
real(8), allocatable :: rev_energy_bins(:)
character(len=MAX_WORD_LEN), allocatable :: names(:)
! Check if mgxs.xml exists
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find mgxs.xml file
call fatal_error("Cross sections XML file '" &
! Could not find MGXS Library file
call fatal_error("Cross sections HDF5 file '" &
// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Reading cross sections XML file...", 5)
call write_message("Reading cross sections HDF5 file...", 5)
! Parse mgxs.xml file
call open_xmldoc(doc, path_cross_sections)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
if (check_for_node(doc, "groups")) then
if (attribute_exists(file_id, "groups")) then
! Get neutron group count
call get_node_value(doc, "groups", energy_groups)
call read_attribute(energy_groups, file_id, "groups")
else
call fatal_error("groups element must exist!")
call fatal_error("'groups' attribute must exist!")
end if
allocate(rev_energy_bins(energy_groups + 1))
allocate(energy_bins(energy_groups + 1))
if (check_for_node(doc, "group_structure")) then
if (attribute_exists(file_id, "group structure")) then
! Get neutron group structure
call get_node_array(doc, "group_structure", energy_bins)
call read_attribute(energy_bins, file_id, "group structure")
else
call fatal_error("group_structures element must exist!")
call fatal_error("'group structure' attribute must exist!")
end if
! First reverse the order of energy_groups
@ -4635,45 +4685,28 @@ contains
energy_bin_avg(i) = HALF * (energy_bins(i) + energy_bins(i + 1))
end do
allocate(inverse_velocities(energy_groups))
if (check_for_node(doc, "inverse_velocities")) then
! Get inverse velocities
call get_node_array(doc, "inverse_velocities", inverse_velocities)
else
! If not given, estimate them by using average energy in group which is
! assumed to be the midpoint
do i = 1, energy_groups
inverse_velocities(i) = ONE / &
(sqrt(TWO * energy_bin_avg(i) / (MASS_NEUTRON_MEV)) * &
C_LIGHT * 100.0_8)
end do
end if
! Get the datasets present in the library
call get_groups(file_id, names)
n_libraries = size(names)
! Get node list of all <xsdata>
call get_node_list(doc, "xsdata", node_xsdata_list)
n_libraries = get_list_size(node_xsdata_list)
! Allocate xs_listings array
! Allocate libraries array
if (n_libraries == 0) then
call fatal_error("At least one <xsdata> element must be present in &
&mgxs.xml file!")
call fatal_error("At least one MGXS data set must be present in &
&mgxs library file!")
else
allocate(libraries(n_libraries))
end if
do i = 1, n_libraries
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i, node_xsdata)
! Get name of material
allocate(libraries(i) % materials(1))
call get_node_value(node_xsdata, "name", libraries(i) % materials(1))
libraries(i) % materials(1) = names(i)
end do
! Close cross sections XML file
call close_xmldoc(doc)
! Close MGXS HDF5 file
call file_close(file_id)
end subroutine read_mg_cross_sections_xml
end subroutine read_mg_cross_sections_header
!===============================================================================
! EXPAND_NATURAL_ELEMENT converts natural elements specified using an <element>
@ -5618,7 +5651,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = ONE
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
end if
! if given weight percent, convert all values so that they are divided
@ -5643,7 +5676,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = ONE
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
end if
x = mat % atom_density(j)
sum_percent = sum_percent + x*awr
@ -5906,67 +5939,6 @@ contains
end do
end subroutine assign_temperatures
!===============================================================================
! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table
! appears at in the model. Later, this list is used to determine the actual
! temperatures to read (which may be different if interpolation is used)
!===============================================================================
subroutine get_temperatures(nuc_temps, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), allocatable, intent(out) :: sab_temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_sab ! index in S(a,b) array
integer :: i_material
real(8) :: temperature ! temperature in Kelvin
allocate(nuc_temps(n_nuclides_total))
allocate(sab_temps(n_sab_tables))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(cells(i) % material(j))
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in nuc_temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
! Add temperature if it hasn't already been added
if (find(nuc_temps(i_nuclide), temperature) == -1) then
call nuc_temps(i_nuclide) % push_back(temperature)
end if
end do NUC_NAMES_LOOP
if (mat % n_sab > 0) then
SAB_NAMES_LOOP: do k = 1, size(mat % sab_names)
! Get index in nuc_temps array
i_sab = sab_dict % get_key(to_lower(mat % sab_names(k)))
! Add temperature if it hasn't already been added
if (find(sab_temps(i_sab), temperature) == -1) then
call sab_temps(i_sab) % push_back(temperature)
end if
end do SAB_NAMES_LOOP
end if
end associate
end do
end do
end subroutine get_temperatures
!===============================================================================
! READ_0K_ELASTIC_SCATTERING
!===============================================================================

139
src/mgxs_data.F90
View file

@ -1,15 +1,17 @@
module mgxs_data
use constants
use constants
use algorithm, only: find
use error, only: fatal_error
use geometry_header, only: get_temperatures
use global
use hdf5_interface
use material_header, only: Material
use mgxs_header
use output, only: write_message
use set_header, only: SetChar
use stl_vector, only: VectorReal
use string, only: to_lower
use xml_interface
implicit none
contains
@ -20,52 +22,41 @@ contains
!===============================================================================
subroutine read_mgxs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_xsdata ! index in <xsdata> list
integer :: i_nuclide ! index in nuclides
character(20) :: name ! name of isotope, e.g. 92235.03c
integer :: i_xsdata ! index in xsdata_dict
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
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
type(Material), pointer :: mat
type(SetChar) :: already_read
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
logical :: get_kfiss, get_fiss
integer :: l
type(DictCharInt) :: xsdata_dict
type(VectorReal), allocatable :: temps(:)
! Check if cross_sections.xml exists
! Check if MGXS Library 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 '" &
! Could not find MGXS Library file
call fatal_error("Cross sections HDF5 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 temperatures
call get_temperatures(cells, materials, material_dict, nuclide_dict, &
n_nuclides_total, temps)
! Get node list of all <xsdata>
call get_node_list(doc, "xsdata", node_xsdata_list)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
! Build dictionary mapping nuclide names to an index in the <xsdata> node
! list
do i = 1, get_list_size(node_xsdata_list)
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i, node_xsdata)
! Get name and create pair (name, i)
call get_node_value(node_xsdata, "name", name)
call xsdata_dict % add_key(to_lower(name), i)
end do
! allocate arrays for ACE table storage and cross section cache
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides_total))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
@ -102,18 +93,22 @@ contains
i_xsdata = xsdata_dict % get_key(to_lower(name))
i_nuclide = mat % nuclide(j)
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i_xsdata, node_xsdata)
call write_message("Loading " // trim(name) // " Data...", 5)
! Check to make sure cross section set exists in the library
if (object_exists(file_id, trim(name))) then
xsdata_group = open_group(file_id, trim(name))
else
call fatal_error("Data for '" // trim(name) // "' does not exist in "&
&// trim(path_cross_sections))
end if
! 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
if (attribute_exists(xsdata_group, "representation")) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
else if (temp_str == 'angle') then
else if (trim(temp_str) == 'angle') then
representation = MGXS_ANGLE
else
call fatal_error("Invalid Data Representation!")
@ -132,8 +127,10 @@ contains
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % init_file(node_xsdata, &
energy_groups, get_kfiss, get_fiss, max_order)
call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
energy_groups, temps(i_nuclide), temperature_method, &
temperature_tolerance, get_kfiss, get_fiss, max_order, &
legendre_to_tabular, legendre_to_tabular_points)
! Add name to dictionary
call already_read % add(name)
@ -172,28 +169,74 @@ contains
subroutine create_macro_xs()
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
integer :: scatt_type
type(VectorReal), allocatable :: kTs(:)
allocate(macro_xs(n_materials))
! Get temperatures to read for each material
call get_mat_kTs(kTs)
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
scatt_type = nuclides_MG(mat % nuclide(1)) % obj % scatt_type
select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
type is (MgxsIso)
allocate(MgxsIso :: macro_xs(i_mat) % obj)
type is (MgxsAngle)
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
call macro_xs(i_mat) % obj % combine(mat, nuclides_MG, energy_groups, &
max_order, scatt_type)
! Do not read materials which we do not actually use in the problem to
! save space
if (allocated(kTs(i_mat) % data)) then
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
energy_groups, max_order, &
temperature_tolerance, &
temperature_method)
end if
end do
end subroutine create_macro_xs
!===============================================================================
! GET_MAT_kTs returns a list of temperatures (in MeV) that each
! material appears at in the model.
!===============================================================================
subroutine get_mat_kTs(kTs)
type(VectorReal), allocatable, intent(out) :: kTs(:)
integer :: i, j
integer :: i_material ! Index in materials array
real(8) :: kT ! temperature in MeV
allocate(kTs(size(materials)))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
kT = cells(i) % sqrtkT(j)**2
else
kT = cells(i) % sqrtkT(1)**2
end if
i_material = material_dict % get_key(cells(i) % material(j))
! Add temperature if it hasn't already been added
if (find(kTs(i_material), kT) == -1) then
call kTs(i_material) % push_back(kT)
end if
end do
end do
end subroutine get_mat_kTs
end module mgxs_data

3131
src/mgxs_header.F90
View file

File diff suppressed because it is too large Load diff

39
src/nuclide_header.F90
View file

@ -3,9 +3,7 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5iget_name_f, h5gget_info_f, &
h5lget_name_by_idx_f, H5_INDEX_NAME_F, H5_ITER_INC_F
use h5lt, only: h5ltpath_valid_f
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: sort, find
use constants
@ -14,7 +12,8 @@ module nuclide_header
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error, only: fatal_error, warning
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape, get_datasets
open_dataset, read_dataset, close_dataset, get_shape, get_datasets, &
object_exists, get_name, get_groups
use list_header, only: ListInt
use math, only: evaluate_legendre
use multipole_header, only: MultipoleArray
@ -196,10 +195,6 @@ module nuclide_header
logical, intent(in) :: master ! if this is the master proc
integer :: i
integer :: storage_type
integer :: max_corder
integer :: n_links
integer :: hdf5_err
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: urr_group, nu_group
@ -210,21 +205,21 @@ module nuclide_header
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(SIZE_T) :: name_len, name_file_len
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp_str
character(MAX_FILE_LEN), allocatable :: dset_names(:)
character(MAX_WORD_LEN), allocatable :: dset_names(:)
character(MAX_WORD_LEN), allocatable :: grp_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
logical :: exists
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
! Get name of nuclide from group
name_len = len(this % name)
call h5iget_name_f(group_id, this % name, name_len, name_file_len, hdf5_err)
this % name = get_name(group_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -336,12 +331,10 @@ module nuclide_header
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
call h5gget_info_f(rxs_group, storage_type, n_links, max_corder, hdf5_err)
do j = 0, n_links - 1
call h5lget_name_by_idx_f(rxs_group, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, &
j, temp_str, hdf5_err, name_len)
if (starts_with(temp_str, "reaction_")) then
call MTs % push_back(int(str_to_int(temp_str(10:12))))
call get_groups(rxs_group, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "reaction_")) then
call MTs % push_back(int(str_to_int(grp_names(j)(10:12))))
end if
end do
@ -357,8 +350,7 @@ module nuclide_header
call close_group(rxs_group)
! Read unresolved resonance probability tables if present
call h5ltpath_valid_f(group_id, 'urr', .true., exists, hdf5_err)
if (exists) then
if (object_exists(group_id, 'urr')) then
this % urr_present = .true.
allocate(this % urr_data(n_temperature))
@ -399,8 +391,7 @@ module nuclide_header
end if
! Check for nu-total
call h5ltpath_valid_f(group_id, 'total_nu', .true., exists, hdf5_err)
if (exists) then
if (object_exists(group_id, 'total_nu')) then
nu_group = open_group(group_id, 'total_nu')
! Read total nu data
@ -419,9 +410,7 @@ module nuclide_header
end if
! Read fission energy release data if present
call h5ltpath_valid_f(group_id, 'fission_energy_release', .true., exists, &
hdf5_err)
if (exists) then
if (object_exists(group_id, 'fission_energy_release')) then
fer_group = open_group(group_id, 'fission_energy_release')
! Check to see if this is polynomial or tabulated data

9
src/output.F90
View file

@ -856,8 +856,13 @@ contains
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
"Total Material"
else
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides(i_nuclide) % name)
if (run_CE) then
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides(i_nuclide) % name)
else
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides_MG(i_nuclide) % obj % name)
end if
end if
indent = indent + 2

12
src/physics.F90
View file

@ -96,10 +96,11 @@ contains
! absorption (including fission)
if (nuc % fissionable) then
call sample_fission(i_nuclide, i_reaction)
if (run_mode == MODE_EIGENVALUE) then
call sample_fission(i_nuclide, i_reaction)
call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
elseif (run_mode == MODE_FIXEDSOURCE) then
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
call sample_fission(i_nuclide, i_reaction)
call create_fission_sites(p, i_nuclide, i_reaction, &
p % secondary_bank, p % n_secondary)
end if
@ -126,6 +127,13 @@ contains
if (.not. p % alive) return
end if
! Kill neutron under certain energy
if (p % E < energy_cutoff) then
p % alive = .false.
p % wgt = ZERO
p % last_wgt = ZERO
end if
end subroutine sample_reaction
!===============================================================================

2
src/physics_mg.F90
View file

@ -73,7 +73,7 @@ contains
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
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
call create_fission_sites(p, p % secondary_bank, p % n_secondary)
end if
end if

17
src/reaction_header.F90
View file

@ -5,7 +5,7 @@ module reaction_header
use constants, only: MAX_WORD_LEN
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape
open_dataset, read_dataset, close_dataset, get_shape, get_groups
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with
@ -42,17 +42,12 @@ contains
integer :: i
integer :: cm
integer :: n_product
integer :: storage_type
integer :: max_corder
integer :: n_links
integer :: hdf5_err
integer(HID_T) :: pgroup
integer(HID_T) :: xs, temp_group
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: j
character(MAX_WORD_LEN) :: name
character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K'
character(MAX_WORD_LEN), allocatable :: grp_names(:)
call read_attribute(this % Q_value, group_id, 'Q_value')
call read_attribute(this % MT, group_id, 'mt')
@ -74,12 +69,10 @@ contains
end do
! Determine number of products
call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err)
n_product = 0
do j = 0, n_links - 1
call h5lget_name_by_idx_f(group_id, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, &
j, name, hdf5_err, name_len)
if (starts_with(name, "product_")) n_product = n_product + 1
call get_groups(group_id, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "product_")) n_product = n_product + 1
end do
! Read products

20
src/sab_header.F90
View file

@ -8,9 +8,9 @@ module sab_header
use distribution_univariate, only: Tabular
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use h5lt, only: h5ltpath_valid_f, h5iget_name_f
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets
open_dataset, read_dataset, close_dataset, get_datasets, object_exists, &
get_name
use secondary_correlated, only: CorrelatedAngleEnergy
use stl_vector, only: VectorInt, VectorReal
use string, only: to_str, str_to_int
@ -66,7 +66,7 @@ module sab_header
end type SabData
type SAlphaBeta
character(100) :: name ! name of table, e.g. lwtr.10t
character(150) :: name ! name of table, e.g. lwtr.10t
real(8) :: awr ! weight of nucleus in neutron masses
real(8), allocatable :: kTs(:) ! temperatures in MeV (k*T)
character(10), allocatable :: nuclides(:) ! List of valid nuclides
@ -92,8 +92,7 @@ contains
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer :: hdf5_err
integer(SIZE_T) :: name_len, name_file_len
integer(SIZE_T) :: name_len
integer(HID_T) :: T_group
integer(HID_T) :: elastic_group
integer(HID_T) :: inelastic_group
@ -103,11 +102,10 @@ contains
integer(HSIZE_T) :: dims3(3)
real(8), allocatable :: temp(:,:)
character(20) :: type
logical :: exists
type(CorrelatedAngleEnergy) :: correlated_dist
character(MAX_WORD_LEN) :: temp_str
character(MAX_FILE_LEN), allocatable :: dset_names(:)
character(MAX_WORD_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
@ -115,7 +113,7 @@ contains
! Get name of table from group
name_len = len(this % name)
call h5iget_name_f(group_id, this % name, name_len, name_file_len, hdf5_err)
this % name = get_name(group_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -207,8 +205,7 @@ contains
T_group = open_group(group_id, temp_str)
! Coherent elastic data
call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err)
if (exists) then
if (object_exists(T_group, 'elastic')) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
@ -250,8 +247,7 @@ contains
end if
! Inelastic data
call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err)
if (exists) then
if (object_exists(T_group, 'inelastic')) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')

416
src/scattdata_header.F90
View file

@ -41,18 +41,20 @@ module scattdata_header
real(8), allocatable :: scattxs(:) ! Isotropic Sigma_{s,g_{in}}
contains
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
procedure(scattdata_init_), deferred :: init ! Initializes ScattData
procedure(scattdata_calc_f_), deferred :: calc_f ! Calculates f, given mu
procedure(scattdata_sample_), deferred :: sample ! sample the scatter event
procedure :: get_matrix => scattdata_get_matrix ! Rebuild scattering matrix
end type ScattData
abstract interface
subroutine scattdata_init_(this, mult, coeffs)
import ScattData
class(ScattData), intent(inout) :: this ! Object to work with
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
subroutine scattdata_init_(this, gmin, gmax, mult, coeffs)
import ScattData, Jagged1D, Jagged2D
class(ScattData), intent(inout) :: this ! Object to work with
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
end subroutine scattdata_init_
pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
@ -79,9 +81,9 @@ module scattdata_header
! Maximal value for rejection sampling from rectangle
type(Jagged1D), allocatable :: max_val(:) ! (Gin % data(Gout))
contains
procedure :: init => scattdatalegendre_init
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
procedure :: init => scattdatalegendre_init
procedure :: calc_f => scattdatalegendre_calc_f
procedure :: sample => scattdatalegendre_sample
end type ScattDataLegendre
type, extends(ScattData) :: ScattDataHistogram
@ -90,10 +92,10 @@ module scattdata_header
! Histogram of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatahistogram_init
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
procedure :: get_matrix => scattdatahistogram_get_matrix
procedure :: init => scattdatahistogram_init
procedure :: calc_f => scattdatahistogram_calc_f
procedure :: sample => scattdatahistogram_sample
procedure :: get_matrix => scattdatahistogram_get_matrix
end type ScattDataHistogram
type, extends(ScattData) :: ScattDataTabular
@ -102,10 +104,10 @@ module scattdata_header
! PDF of f(mu) (dist has CDF)
type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
contains
procedure :: init => scattdatatabular_init
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
procedure :: get_matrix => scattdatatabular_get_matrix
procedure :: init => scattdatatabular_init
procedure :: calc_f => scattdatatabular_calc_f
procedure :: sample => scattdatatabular_sample
procedure :: get_matrix => scattdatatabular_get_matrix
end type ScattDataTabular
!===============================================================================
@ -122,13 +124,15 @@ 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(inout) :: energy(:, :) ! Energy Transfer Matrix
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
subroutine scattdata_init(this, order, gmin, gmax, energy, mult)
class(ScattData), intent(inout) :: this ! Object to work on
integer, intent(in) :: order ! Data Order
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(inout) :: energy(:) ! Energy Transfer Matrix
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
integer :: groups, gmin, gmax, gin
integer :: groups, gin
real(8) :: norm
groups = size(energy, dim=1)
@ -138,83 +142,80 @@ contains
allocate(this % energy(groups))
allocate(this % mult(groups))
allocate(this % dist(groups))
! Use energy to find the gmin and gmax values
! Also set energy values when doing it
this % gmin = gmin
this % gmax = gmax
! Set the outgoing energy PDF values
do gin = 1, groups
! Make sure energy is normalized (i.e., CDF is 1)
norm = sum(energy(:, gin))
if (norm /= ZERO) energy(:, gin) = energy(:, gin) / norm
! Find gmin by checking the P0 moment
do gmin = 1, groups
if (energy(gmin, gin) > ZERO) exit
end do
! Find gmax by checking the P0 moment
do gmax = groups, 1, -1
if (energy(gmax, gin) > ZERO) exit
end do
! Treat the case of all zeros
if (gmin > gmax) then
gmin = gin
gmax = gin
! By not changing energy(gin) here we are leaving it as zero
end if
allocate(this % energy(gin) % data(gmin:gmax))
this % energy(gin) % data(gmin:gmax) = energy(gmin:gmax, gin)
allocate(this % mult(gin) % data(gmin:gmax))
this % mult(gin) % data(gmin:gmax) = mult(gmin:gmax, gin)
allocate(this % dist(gin) % data(order, gmin:gmax))
norm = sum(energy(gin) % data(:))
if (norm /= ZERO) energy(gin) % data(:) = energy(gin) % data(:) / norm
! Set the values
allocate(this % energy(gin) % data(gmin(gin):gmax(gin)))
this % energy(gin) % data(:) = energy(gin) % data(:)
allocate(this % mult(gin) % data(gmin(gin):gmax(gin)))
this % mult(gin) % data(gmin(gin):gmax(gin)) = &
mult(gin) % data(gmin(gin):gmax(gin))
allocate(this % dist(gin) % data(order, gmin(gin):gmax(gin)))
this % dist(gin) % data = ZERO
this % gmin(gin) = gmin
this % gmax(gin) = gmax
end do
end subroutine scattdata_init
subroutine scattdatalegendre_init(this, mult, coeffs)
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
subroutine scattdatalegendre_init(this, gmin, gmax, mult, coeffs)
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
real(8) :: dmu, mu, f, norm
integer :: imu, Nmu, gout, gin, groups, order
real(8), allocatable :: energy(:, :)
real(8), allocatable :: matrix(:, :, :)
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs, dim=3)
order = size(coeffs, dim=1)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(order, groups, groups))
matrix (:, :, :)= coeffs
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data = coeffs(gin) % data
end do
! Get scattxs value
allocate(this % scattxs(groups))
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
this % scattxs(:) = sum(matrix(1, :, :), dim=1)
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
end do
allocate(energy(groups, groups))
energy(:, :) = ZERO
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
do gout = 1, groups
norm = matrix(1, gout, gin)
energy(gout, gin) = norm
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
energy(gin) % data = ZERO
do gout = gmin(gin), gmax(gin)
norm = matrix(gin) % data(1, gout)
energy(gin) % data(gout) = norm
if (norm /= ZERO) then
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
end if
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init(this, order, gmin, gmax, energy, mult)
allocate(this % max_val(groups))
! Set dist values from matrix and initialize max_val
do gin = 1, groups
do gout = this % gmin(gin), this % gmax(gin)
this % dist(gin) % data(:, gout) = matrix(:, gout, gin)
do gout = gmin(gin), gmax(gin)
this % dist(gin) % data(:, gout) = matrix(gin) % data(:, gout)
end do
allocate(this % max_val(gin) % data(this % gmin(gin):this % gmax(gin)))
allocate(this % max_val(gin) % data(gmin(gin):gmax(gin)))
this % max_val(gin) % data(:) = ZERO
end do
@ -223,7 +224,7 @@ contains
Nmu = 1001
dmu = TWO / real(Nmu - 1, 8)
do gin = 1, groups
do gout = this % gmin(gin), this % gmax(gin)
do gout = gmin(gin), gmax(gin)
do imu = 1, Nmu
! Update mu. Do first and last seperate to avoid float errors
if (imu == 1) then
@ -246,44 +247,51 @@ contains
end do
end subroutine scattdatalegendre_init
subroutine scattdatahistogram_init(this, mult, coeffs)
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
subroutine scattdatahistogram_init(this, gmin, gmax, mult, coeffs)
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
integer :: imu, gin, gout, groups, order
real(8) :: norm
real(8), allocatable :: energy(:, :)
real(8), allocatable :: matrix(:, :, :)
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs, dim=3)
order = size(coeffs, dim=1)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(order, groups, groups))
matrix(:, :, :) = coeffs
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data = coeffs(gin) % data
end do
! Get scattxs value
allocate(this % scattxs(groups))
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
this % scattxs(:) = sum(sum(matrix(:, :, :), dim=1), dim=1)
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
end do
allocate(energy(groups, groups))
energy(:, :) = ZERO
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = 1, groups
norm = sum(matrix(:, gout, gin))
energy(gout, gin) = norm
norm = sum(matrix(gin) % data(:, gout))
energy(gin) % data(gout) = norm
if (norm /= ZERO) then
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
end if
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init(this, order, gmin, gmax, energy, mult)
allocate(this % mu(order))
this % dmu = TWO / real(order, 8)
@ -296,17 +304,16 @@ contains
! also saving the original histogram in fmu
allocate(this % fmu(groups))
do gin = 1, groups
allocate(this % fmu(gin) % data(order, &
this % gmin(gin):this % gmax(gin)))
do gout = this % gmin(gin), this % gmax(gin)
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
! Store the histogram
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
! Integrate the histogram
this % dist(gin) % data(1, gout) = &
this % dmu * matrix(1, gout, gin)
this % dmu * matrix(gin) % data(1, gout)
do imu = 2, order
this % dist(gin) % data(imu, gout) = &
this % dmu * matrix(imu, gout, gin) + &
this % dmu * matrix(gin) % data(imu, gout) + &
this % dist(gin) % data(imu - 1, gout)
end do
@ -323,22 +330,27 @@ contains
end subroutine scattdatahistogram_init
subroutine scattdatatabular_init(this, mult, coeffs)
class(ScattDataTabular), intent(inout) :: this ! Object to work on
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
subroutine scattdatatabular_init(this, gmin, gmax, mult, coeffs)
class(ScattDataTabular), intent(inout) :: this ! Object to work on
integer, intent(in) :: gmin(:) ! Min Gout
integer, intent(in) :: gmax(:) ! Max Gout
type(Jagged1D), intent(in) :: mult(:) ! Scatter Prod'n Matrix
type(Jagged2D), intent(in) :: coeffs(:) ! Coefficients to use
integer :: imu, gin, gout, groups, order
real(8) :: norm
real(8), allocatable :: energy(:, :)
real(8), allocatable :: matrix(:, :, :)
type(Jagged1D), allocatable :: energy(:)
type(Jagged2D), allocatable :: matrix(:)
groups = size(coeffs, dim=3)
order = size(coeffs, dim=1)
groups = size(coeffs)
order = size(coeffs(1) % data, dim=1)
! make a copy of coeffs that we can use to extract data and normalize
allocate(matrix(order, groups, groups))
matrix(:, :, :) = coeffs
allocate(matrix(groups))
do gin = 1, groups
allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
matrix(gin) % data = coeffs(gin) % data
end do
! Build the angular distribution mu values
allocate(this % mu(order))
@ -356,39 +368,40 @@ contains
! over all outgoing groups
do gin = 1, groups
norm = ZERO
do gout = 1, groups
do gout = gmin(gin), gmax(gin)
do imu = 2, order
norm = norm + HALF * this % dmu * (matrix(imu - 1, gout, gin) + &
matrix(imu, gout, gin))
norm = norm + HALF * this % dmu * &
(matrix(gin) % data(imu - 1, gout) + &
matrix(gin) % data(imu, gout))
end do
end do
this % scattxs(gin) = norm
end do
allocate(energy(groups, groups))
energy(:, :) = ZERO
allocate(energy(groups))
! Build energy transfer probability matrix from data in matrix
do gin = 1, groups
do gout = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
norm = ZERO
do imu = 2, order
norm = norm + HALF * this % dmu * &
(matrix(imu - 1, gout, gin) + matrix(imu, gout, gin))
(matrix(gin) % data(imu - 1, gout) + &
matrix(gin) % data(imu, gout))
end do
energy(gout, gin) = norm
energy(gin) % data(gout) = norm
end do
end do
call scattdata_init(this, order, energy, mult)
call scattdata_init(this, order, gmin, gmax, energy, mult)
! Calculate f(mu) and integrate it so we can avoid rejection sampling
allocate(this % fmu(groups))
do gin = 1, groups
allocate(this % fmu(gin) % data(order, &
this % gmin(gin):this % gmax(gin)))
do gout = this % gmin(gin), this % gmax(gin)
allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
! Coeffs contain f(mu), put in f(mu) as that is where the
! PDF lives
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
! Force positivity
do imu = 1, order
@ -640,10 +653,10 @@ contains
! using ScattData's information of fmu/dist, energy, and scattxs
!===============================================================================
pure function scattdata_get_matrix(this, req_order) result(matrix)
class(ScattData), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
subroutine scattdata_get_matrix(this, req_order, matrix)
class(ScattData), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
@ -651,66 +664,189 @@ contains
! Set gin and gout for getting the order
order = min(req_order, size(this % dist(1) % data, dim=1))
allocate(matrix(order, groups, groups))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
matrix(:, :, :) = ZERO
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(:, gout, gin) = this % scattxs(gin) * &
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % dist(gin) % data(1:order, gout)
end do
end do
end function scattdata_get_matrix
end subroutine scattdata_get_matrix
pure function scattdatahistogram_get_matrix(this, req_order) result(matrix)
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
subroutine scattdatahistogram_get_matrix(this, req_order, matrix)
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
groups = size(this % energy)
order = min(req_order, size(this % dist(1) % data, dim=1))
allocate(matrix(order, groups, groups))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
matrix(:, :, :) = ZERO
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(:, gout, gin) = this % scattxs(gin) * &
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % fmu(gin) % data(1:order, gout)
end do
end do
end function scattdatahistogram_get_matrix
end subroutine scattdatahistogram_get_matrix
pure function scattdatatabular_get_matrix(this, req_order) result(matrix)
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
subroutine scattdatatabular_get_matrix(this, req_order, matrix)
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
integer, intent(in) :: req_order ! Requested order of matrix
type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
integer :: order, groups, gin, gout
groups = size(this % energy)
order = min(req_order, size(this % dist(1) % data, dim=1))
allocate(matrix(order, groups, groups))
if (allocated(matrix)) deallocate(matrix)
allocate(matrix(groups))
! Initialize to 0; this way the zero entries in the dense matrix dont
! need to be explicitly set, requiring a significant increase in the
! lines of code.
matrix(:, :, :) = ZERO
do gin = 1, groups
allocate(matrix(gin) % data(order, groups))
do gout = this % gmin(gin), this % gmax(gin)
matrix(:, gout, gin) = this % scattxs(gin) * &
matrix(gin) % data(:, gout) = this % scattxs(gin) * &
this % energy(gin) % data(gout) * &
this % fmu(gin) % data(1:order, gout)
end do
end do
end function scattdatatabular_get_matrix
end subroutine scattdatatabular_get_matrix
!===============================================================================
! JAGGED_FROM_DENSE_*D Creates a jagged array from a sparse dense matrix.
! The user can supply a key which indicates the values to remove, but the
! default is ZERO
!===============================================================================
subroutine jagged_from_dense_1D(dense, jagged, lo_bounds_, hi_bounds_, key_)
real(8), intent(in) :: dense(:, :)
type(Jagged1D), allocatable, intent(inout) :: jagged(:)
real(8), intent(in), optional :: key_
integer, intent(inout), allocatable, optional :: lo_bounds_(:)
integer, intent(inout), allocatable, optional :: hi_bounds_(:)
real(8) :: key
integer :: i, jmin, jmax
integer, allocatable :: lo_bounds(:), hi_bounds(:)
if (present(key_)) then
key = key_
else
key = ZERO
end if
allocate(lo_bounds(size(dense, dim=2)))
allocate(hi_bounds(size(dense, dim=2)))
if (allocated(jagged)) deallocate(jagged)
allocate(jagged(size(dense, dim=2)))
do i = 1, size(dense, dim=2)
! Find the min and max j values
do jmin = 1, size(dense, dim=1)
if (dense(jmin, i) /= key) exit
end do
do jmax = size(dense, dim=1), 1, -1
if (dense(jmax, i) /= key) exit
end do
! Treat the case of all values matching the key
if (jmin > jmax) then
jmin = i
jmax = i
end if
! Now store the jagged row
allocate(jagged(i) % data(jmin:jmax))
jagged(i) % data(jmin:jmax) = dense(jmin:jmax, i)
lo_bounds(i) = jmin
hi_bounds(i) = jmax
end do
if (present(lo_bounds_)) then
if (allocated(lo_bounds_)) deallocate(lo_bounds_)
allocate(lo_bounds_(size(dense, dim=2)))
lo_bounds_ = lo_bounds
end if
if (present(hi_bounds_)) then
if (allocated(hi_bounds_)) deallocate(hi_bounds_)
allocate(hi_bounds_(size(dense, dim=2)))
hi_bounds_ = hi_bounds
end if
end subroutine jagged_from_dense_1D
subroutine jagged_from_dense_2D(dense, jagged, lo_bounds_, hi_bounds_, key_)
real(8), intent(in) :: dense(:, :, :)
type(Jagged2D), allocatable, intent(inout) :: jagged(:)
real(8), intent(in), optional :: key_
integer, intent(inout), allocatable, optional :: lo_bounds_(:)
integer, intent(inout), allocatable, optional :: hi_bounds_(:)
real(8) :: key
integer :: i, jmin, jmax
integer, allocatable :: lo_bounds(:), hi_bounds(:)
if (present(key_)) then
key = key_
else
key = ZERO
end if
allocate(lo_bounds(size(dense, dim=3)))
allocate(hi_bounds(size(dense, dim=3)))
if (allocated(jagged)) deallocate(jagged)
allocate(jagged(size(dense, dim=3)))
do i = 1, size(dense, dim=3)
! Find the min and max j values
do jmin = 1, size(dense, dim=2)
if (any(dense(:, jmin, i) /= key)) exit
end do
do jmax = size(dense, dim=2), 1, -1
if (any(dense(:, jmax, i) /= key)) exit
end do
! Treat the case of all values matching the key
if (jmin > jmax) then
jmin = i
jmax = i
end if
! Now store the jagged row
allocate(jagged(i) % data(size(dense, dim=1), jmin:jmax))
jagged(i) % data(:, jmin:jmax) = dense(:, jmin:jmax, i)
lo_bounds(i) = jmin
hi_bounds(i) = jmax
end do
if (present(lo_bounds_)) then
if (allocated(lo_bounds_)) deallocate(lo_bounds_)
allocate(lo_bounds_(size(dense, dim=3)))
lo_bounds_ = lo_bounds
end if
if (present(hi_bounds_)) then
if (allocated(hi_bounds_)) deallocate(hi_bounds_)
allocate(hi_bounds_(size(dense, dim=3)))
hi_bounds_ = hi_bounds
end if
end subroutine jagged_from_dense_2D
end module scattdata_header

18
src/secondary_uncorrelated.F90
View file

@ -1,6 +1,5 @@
module secondary_uncorrelated
use h5lt, only: h5ltpath_valid_f
use hdf5, only: HID_T
use angle_distribution, only: AngleDistribution
@ -9,7 +8,8 @@ module secondary_uncorrelated
use energy_distribution, only: EnergyDistribution, LevelInelastic, &
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton
use error, only: warning
use hdf5_interface, only: read_attribute, open_group, close_group
use hdf5_interface, only: read_attribute, open_group, close_group, &
object_exists
use random_lcg, only: prn
!===============================================================================
@ -56,25 +56,19 @@ contains
class(UncorrelatedAngleEnergy), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
logical :: exists
integer :: hdf5_err
integer(HID_T) :: energy_group
integer(HID_T) :: angle_group
character(MAX_WORD_LEN) :: type
! Check if energy group is present
call h5ltpath_valid_f(group_id, 'angle', .true., exists, hdf5_err)
if (exists) then
! Check if angle group is present & read
if (object_exists(group_id, 'angle')) then
angle_group = open_group(group_id, 'angle')
call this%angle%from_hdf5(angle_group)
call close_group(angle_group)
end if
! Check if energy group is present
call h5ltpath_valid_f(group_id, 'energy', .true., exists, hdf5_err)
if (exists) then
! Check if energy group is present & read
if (object_exists(group_id, 'energy')) then
energy_group = open_group(group_id, 'energy')
call read_attribute(type, energy_group, 'type')
select case (type)

17
src/state_point.F90
View file

@ -251,11 +251,20 @@ contains
allocate(str_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (tally % nuclide_bins(j) > 0) then
i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
if (run_CE) then
i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
else
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
else
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
i_xs = index(nuclides_MG(tally % nuclide_bins(j)) % obj % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name(1 : i_xs-1)
else
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name
end if
end if
else
str_array(j) = 'total'

43
src/tally.F90
View file

@ -1149,6 +1149,8 @@ contains
! Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide > 0) then
nucxs => nuclides_MG(i_nuclide) % obj
! And since we haven't calculated this temperature index yet, do so now
call nucxs % find_temperature(p % sqrtkT)
end if
i = 0
@ -1228,11 +1230,21 @@ contains
else
score = p % last_wgt
end if
score = score * inverse_velocities(p_g) / material_xs % total * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) / &
matxs % get_xs('total', p_g, UVW=p_uvw) * flux
else
score = matxs % get_xs('inv_vel', p_g, UVW=p_uvw) * flux
end if
else
! For inverse velocity, we need no cross section
score = flux * inverse_velocities(p_g)
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inv_vel', p_g, UVW=p_uvw) * &
atom_density * flux
else
score = flux * matxs % get_xs('inv_vel', p_g, UVW=p_uvw)
end if
end if
@ -1853,17 +1865,24 @@ contains
i_nuclide = t % nuclide_bins(k)
! Check to see if this nuclide was in the material of our collision.
do m = 1, mat % n_nuclides
if (mat % nuclide(m) == i_nuclide) then
atom_density = mat % atom_density(m)
exit
end if
end do
if (i_nuclide > 0) then
atom_density = -ONE
! Check to see if this nuclide was in the material of our collision
do m = 1, mat % n_nuclides
if (mat % nuclide(m) == i_nuclide) then
atom_density = mat % atom_density(m)
exit
end if
end do
else
atom_density = ZERO
end if
! Determine score for each bin
call score_general(p, t, (k-1)*t % n_score_bins, filter_index, &
i_nuclide, atom_density, filter_weight)
! If we found the nuclide, determine the score for each bin
if (atom_density >= ZERO) then
call score_general(p, t, (k-1)*t % n_score_bins, filter_index, &
i_nuclide, atom_density, filter_weight)
end if
end do NUCLIDE_LOOP

3
src/tracking.F90
View file

@ -92,6 +92,9 @@ contains
! 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
! Update the temperature index
call macro_xs(p % material) % obj % find_temperature(p % sqrtkT)
! Get the data
call macro_xs(p % material) % obj % calculate_xs(p % g, &
p % coord(p % n_coord) % uvw, material_xs)
else

BIN
tests/1d_mgxs.h5 Normal file
View file

Binary file not shown.

9859
tests/1d_mgxs.xml
View file

File diff suppressed because it is too large Load diff

92
tests/input_set.py
View file

@ -1,8 +1,9 @@
import numpy as np
import openmc
from openmc.source import Source
from openmc.stats import Box
import numpy as np
class InputSet(object):
def __init__(self):
@ -728,60 +729,59 @@ class AssemblyInputSet(object):
class MGInputSet(InputSet):
def build_default_materials_and_geometry(self):
def build_default_materials_and_geometry(self, reps=None, as_macro=True):
# Define materials needed for 1D/1G slab problem
uo2_data = openmc.Macroscopic('uo2_iso')
uo2 = openmc.Material(name='UO2', material_id=1)
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
mat_names = ['uo2', 'clad', 'lwtr']
mgxs_reps = ['ang', 'ang_mu', 'iso', 'iso_mu']
clad_data = openmc.Macroscopic('clad_ang_mu')
clad = openmc.Material(name='Clad', material_id=2)
clad.set_density('macro', 1.0)
clad.add_macroscopic(clad_data)
if reps is None:
reps = mgxs_reps
water_data = openmc.Macroscopic('lwtr_iso_mu')
water = openmc.Material(name='LWTR', material_id=3)
water.set_density('macro', 1.0)
water.add_macroscopic(water_data)
xs = []
mats = []
i = 0
for mat in mat_names:
for rep in reps:
i += 1
if as_macro:
xs.append(openmc.Macroscopic(mat + '_' + rep))
mats.append(openmc.Material(name=str(i)))
mats[-1].set_density('macro', 1.)
mats[-1].add_macroscopic(xs[-1])
else:
xs.append(openmc.Nuclide(mat + '_' + rep))
mats.append(openmc.Material(name=str(i)))
mats[-1].set_density('atom/b-cm', 1.)
mats[-1].add_nuclide(xs[-1].name, 1.0, 'ao')
# Define the materials file.
self.materials += (uo2, clad, water)
# Define the materials file
self.xs_data = xs
self.materials += mats
# 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')
left = openmc.XPlane(x0=0.0, boundary_type='reflective')
right = openmc.XPlane(x0=10.0, boundary_type='reflective')
bottom = openmc.YPlane(y0=0.0, boundary_type='reflective')
top = openmc.YPlane(y0=10.0, boundary_type='reflective')
# for each material add a plane
planes = [openmc.ZPlane(z0=0.0, boundary_type='reflective')]
dz = round(5. / float(len(mats)), 4)
for i in range(len(mats) - 1):
planes.append(openmc.ZPlane(z0=dz * float(i + 1)))
planes.append(openmc.ZPlane(z0=5.0, 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 cells for each material
cells = []
xy = +left & -right & +bottom & -top
for i, mat in enumerate(mats):
cells.append(openmc.Cell())
cells[-1].region = xy & +planes[i] & -planes[i + 1]
cells[-1].fill = mat
# Define root universe.
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cells((c1, c2, c3))
root.add_cells(cells)
# Assign root universe to geometry
self.geometry.root_universe = root
@ -791,9 +791,9 @@ class MGInputSet(InputSet):
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]))
[10.0, 10.0, 5.]))
self.settings.energy_mode = "multi-group"
self.settings.cross_sections = "../1d_mgxs.xml"
self.settings.cross_sections = "../1d_mgxs.h5"
def build_defualt_plots(self):
plot = openmc.Plot()

2
tests/run_tests.py
View file

@ -117,7 +117,7 @@ def cleanup(path):
for dirpath, dirnames, filenames in os.walk(path):
for fname in filenames:
for ext in ['.h5', '.ppm', '.voxel']:
if fname.endswith(ext):
if fname.endswith(ext) and fname != '1d_mgxs.h5':
os.remove(os.path.join(dirpath, fname))

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

@ -0,0 +1 @@
c3581501d1486293c255390251d20584431a198fbb7c07dc2879dc95dc96e26786d3913ce0db8007f542468fd137ff3267824844c03b4687fdad81ea568c92d7

3
tests/test_create_fission_neutrons/results_true.dat Normal file
View file

@ -0,0 +1,3 @@
tally 1:
sum = 2.056839E+02
sum_sq = 4.244628E+03

80
tests/test_create_fission_neutrons/test_create_fission_neutrons.py Executable file
View file

@ -0,0 +1,80 @@
#!/usr/bin/env python
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
import openmc
class CreateFissionNeutronsTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Material is composed of H-1 and U-235
mat = openmc.Material(material_id=1, name='mat')
mat.set_density('atom/b-cm', 0.069335)
mat.add_nuclide('H1', 40.0)
mat.add_nuclide('U235', 1.0)
materials_file = openmc.Materials([mat])
materials_file.export_to_xml()
# Cell is box with reflective boundary
x1 = openmc.XPlane(surface_id=1, x0=-1)
x2 = openmc.XPlane(surface_id=2, x0=1)
y1 = openmc.YPlane(surface_id=3, y0=-1)
y2 = openmc.YPlane(surface_id=4, y0=1)
z1 = openmc.ZPlane(surface_id=5, z0=-1)
z2 = openmc.ZPlane(surface_id=6, z0=1)
for surface in [x1, x2, y1, y2, z1, z2]:
surface.boundary_type = 'reflective'
box = openmc.Cell(cell_id=1, name='box')
box.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2
box.fill = mat
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cell(box)
geometry = openmc.Geometry(root)
geometry.export_to_xml()
# Set the running parameters
settings_file = openmc.Settings()
settings_file.run_mode = 'fixed source'
settings_file.batches = 10
settings_file.particles = 100
settings_file.create_fission_neutrons = False
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
watt_dist = openmc.stats.Watt()
settings_file.source = openmc.source.Source(space=uniform_dist,
energy=watt_dist)
settings_file.export_to_xml()
# Create tallies
tallies = openmc.Tallies()
tally = openmc.Tally(1)
tally.scores = ['flux']
tallies.append(tally)
tallies.export_to_xml()
def _get_results(self):
"""Digest info in the statepoint and return as a string."""
# Read the statepoint file.
sp = openmc.StatePoint(self._sp_name)
# Write out tally data.
outstr = ''
t = sp.get_tally()
outstr += 'tally {0}:\n'.format(t.id)
outstr += 'sum = {0:12.6E}\n'.format(t.sum[0, 0, 0])
outstr += 'sum_sq = {0:12.6E}\n'.format(t.sum_sq[0, 0, 0])
return outstr
def _cleanup(self):
super(CreateFissionNeutronsTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
harness = CreateFissionNeutronsTestHarness('statepoint.10.h5', True)
harness.main()

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

@ -0,0 +1 @@
4dc6a7b131f6757ecc9b06e93d425712f0813c546ed9eaa0aafa4990383b2f3a5b742a5a39ef5f27300d2c861b344e78de2d2c569b9439e7893dc06b03bf3b02

3
tests/test_energy_cutoff/results_true.dat Normal file
View file

@ -0,0 +1,3 @@
tally 1:
sum = 0.000000E+00
sum_sq = 0.000000E+00

84
tests/test_energy_cutoff/test_energy_cutoff.py Executable 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 PyAPITestHarness
import openmc
class EnergyCutoffTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Set energy cutoff
energy_cutoff = 4e-6
# Material is composed of H-1
mat = openmc.Material(material_id=1, name='mat')
mat.set_density('atom/b-cm', 0.069335)
mat.add_nuclide('H1', 40.0)
materials_file = openmc.Materials([mat])
materials_file.export_to_xml()
# Cell is box with reflective boundary
x1 = openmc.XPlane(surface_id=1, x0=-1)
x2 = openmc.XPlane(surface_id=2, x0=1)
y1 = openmc.YPlane(surface_id=3, y0=-1)
y2 = openmc.YPlane(surface_id=4, y0=1)
z1 = openmc.ZPlane(surface_id=5, z0=-1)
z2 = openmc.ZPlane(surface_id=6, z0=1)
for surface in [x1, x2, y1, y2, z1, z2]:
surface.boundary_type = 'reflective'
box = openmc.Cell(cell_id=1, name='box')
box.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2
box.fill = mat
root = openmc.Universe(universe_id=0, name='root universe')
root.add_cell(box)
geometry = openmc.Geometry(root)
geometry.export_to_xml()
# Set the running parameters
settings_file = openmc.Settings()
settings_file.run_mode = 'fixed source'
settings_file.batches = 10
settings_file.particles = 100
settings_file.cutoff = {'energy': energy_cutoff}
bounds = [-1, -1, -1, 1, 1, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
watt_dist = openmc.stats.Watt()
settings_file.source = openmc.source.Source(space=uniform_dist,
energy=watt_dist)
settings_file.export_to_xml()
# Tally flux under energy cutoff
tallies = openmc.Tallies()
tally = openmc.Tally(1)
tally.scores = ['flux']
energy_filter = openmc.filter.EnergyFilter((0.0, energy_cutoff))
tally.filters = [energy_filter]
tallies.append(tally)
tallies.export_to_xml()
def _get_results(self):
"""Digest info in the statepoint and return as a string."""
# Read the statepoint file.
sp = openmc.StatePoint(self._sp_name)
# Write out tally data.
outstr = ''
t = sp.get_tally()
outstr += 'tally {0}:\n'.format(t.id)
outstr += 'sum = {0:12.6E}\n'.format(t.sum[0, 0, 0])
outstr += 'sum_sq = {0:12.6E}\n'.format(t.sum_sq[0, 0, 0])
return outstr
def _cleanup(self):
super(EnergyCutoffTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'tallies.xml')
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
harness = EnergyCutoffTestHarness('statepoint.10.h5', True)
harness.main()

2
tests/test_mg_basic/inputs_true.dat
View file

@ -1 +1 @@
e843dbee8b989142d68e78f7a3e83309a9982c0127974ba4b6a58c823ce298e8fffc8828b48aa696f045818a2a80ef4003f5a191f33aacef7d7dd72c950855a6
115218221500e60ea43145875324a9a6800366fcc97357ef4cd8182f9253e114cbe1c7481dd3bdd8d7dce5983e2cb757ae0a2e05ede236d94cba7909703ba7f7

2
tests/test_mg_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.047136E+00 2.765964E-02
1.160653E+00 6.227102E-02

2
tests/test_mg_max_order/inputs_true.dat
View file

@ -1 +1 @@
7d508b1f3a2661566b8e8cb76fee61aecb96e8b60d633b06f13c4600bd854ea3366cdebd033c0a71ffb8adb90a9aeb64fe5ac0ef3235260921f5689c93e54305
139d760cd83eab001ed3fd52d7146fb5f30b7021b26ada6a1e425f5446185ac6fa61ac2ac58aa9e895981e10540160b7ae7428c833024679ef1ca762715e2d51

2
tests/test_mg_max_order/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.102093E+00 4.962190E-02
1.098323E+00 1.979833E-02

75
tests/test_mg_max_order/test_mg_max_order.py
View file

@ -2,83 +2,26 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
from testing_harness import 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')
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')
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')
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 += (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
class MGMaxOrderTestHarness(PyAPITestHarness):
def __init__(self, statepoint_name, tallies_present, mg=False):
PyAPITestHarness.__init__(self, statepoint_name, tallies_present)
self._input_set = MGNuclideInputSet()
self._input_set = MGInputSet()
def _build_inputs(self):
"""Write input XML files."""
reps = ['iso']
self._input_set.build_default_materials_and_geometry(reps=reps)
self._input_set.build_default_settings()
# Set P1 scattering
self._input_set.settings.max_order = 1
# Call standard input build
super(MGMaxOrderTestHarness, self)._build_inputs()
self._input_set.export()
if __name__ == '__main__':
harness = MGMaxOrderTestHarness('statepoint.10.*', False, mg=True)

2
tests/test_mg_nuclide/inputs_true.dat
View file

@ -1 +1 @@
be296da93031694b2915e1a10e3e6fd663d612cdd29a84745e3ccb0065c088b4bcda45f6919962f6ee784efe52ebe178bf7d3ce019859637ae57c2da1e240a04
6d03b988671fe21dfe2a908685536edb3059e5bdda7974ddfa65c7f63cecb7f6b3fce49bc1533ef2787e284121779653cb14a66cfd709258e4b07d8bd7571843

2
tests/test_mg_nuclide/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.047136E+00 2.765964E-02
1.160653E+00 6.227102E-02

73
tests/test_mg_nuclide/test_mg_nuclide.py
View file

@ -2,81 +2,22 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
from testing_harness import 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('sum', 1.0)
uo2.add_nuclide("uo2_iso", 1.0)
clad = openmc.Material(name='Clad', material_id=2)
clad.set_density('sum', 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('sum', 1.0)
water.add_nuclide("lwtr_iso_mu", 1.0)
# Define the materials file.
self.materials.default_xs = '71c'
self.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
class MGNuclideTestHarness(PyAPITestHarness):
def __init__(self, statepoint_name, tallies_present, mg=False):
PyAPITestHarness.__init__(self, statepoint_name, tallies_present)
self._input_set = MGNuclideInputSet()
self._input_set = MGInputSet()
def _build_inputs(self):
super(MGNuclideTestHarness, self)._build_inputs()
"""Write input XML files."""
self._input_set.build_default_materials_and_geometry(as_macro=False)
self._input_set.build_default_settings()
self._input_set.export()
if __name__ == '__main__':

2
tests/test_mg_tallies/inputs_true.dat
View file

@ -1 +1 @@
b607875dcaecb7110e396a62100182818b8b2853ec9921194b7ab00d6156373b01394ad8bd729babf7d258e6d9c599f6944c6857cd9d2f65ea98e245ee3cf010
8c496e13f9456bedcc0789ffdca81913596cca3c7867717ced38c3063b6cdc8949bd204c98723b72c1b456e3fc85fc90e6710482b0c74e5e29283902a5d9ed50

2907
tests/test_mg_tallies/results_true.dat
View file

File diff suppressed because it is too large Load diff

74
tests/test_mg_tallies/test_mg_tallies.py
View file

@ -3,38 +3,84 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness, PyAPITestHarness
from testing_harness import HashedPyAPITestHarness
import openmc
class MGTalliesTestHarness(PyAPITestHarness):
class MGTalliesTestHarness(HashedPyAPITestHarness):
def _build_inputs(self):
"""Write input XML files."""
self._input_set.build_default_materials_and_geometry(as_macro=False)
self._input_set.build_default_settings()
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [17, 17, 1]
mesh.dimension = [1, 1, 10]
mesh.lower_left = [0.0, 0.0, 0.0]
mesh.upper_right = [21.42, 21.42, 100.0]
mesh.upper_right = [10, 10, 5]
# Instantiate some tally filters
energy_filter = openmc.EnergyFilter([0.0, 20.0])
energyout_filter = openmc.EnergyoutFilter([0.0, 20.0])
mesh_filter = openmc.MeshFilter(mesh)
mat_filter = openmc.MaterialFilter([1,2,3])
mat_ids = [mat.id for mat in self._input_set.materials]
mat_filter = openmc.MaterialFilter(mat_ids)
tally1 = openmc.Tally(tally_id=1)
tally1.filters = [mesh_filter]
tally1.scores = ['total', 'absorption', 'flux',
'fission', 'nu-fission']
nuclides = [xs.name for xs in self._input_set.xs_data]
tally2 = openmc.Tally(tally_id=2)
tally2.filters = [mat_filter, energy_filter, energyout_filter]
tally2.scores = ['scatter', 'nu-scatter']
scores= {False: ['total', 'absorption', 'flux', 'fission', 'nu-fission'],
True: ['total', 'absorption', 'fission', 'nu-fission']}
self._input_set.tallies = openmc.Tallies([tally1, tally2])
tallies = []
for do_nuclides in [False, True]:
tallies.append(openmc.Tally())
tallies[-1].filters = [mesh_filter]
tallies[-1].estimator = 'analog'
tallies[-1].scores = scores[do_nuclides]
if do_nuclides:
tallies[-1].nuclides = nuclides
super(MGTalliesTestHarness, self)._build_inputs()
tallies.append(openmc.Tally())
tallies[-1].filters = [mesh_filter]
tallies[-1].estimator = 'tracklength'
tallies[-1].scores = scores[do_nuclides]
if do_nuclides:
tallies[-1].nuclides = nuclides
tallies.append(openmc.Tally())
tallies[-1].filters = [mat_filter, energy_filter]
tallies[-1].estimator = 'analog'
tallies[-1].scores = scores[do_nuclides] + ['scatter',
'nu-scatter']
if do_nuclides:
tallies[-1].nuclides = nuclides
tallies.append(openmc.Tally())
tallies[-1].filters = [mat_filter, energy_filter]
tallies[-1].estimator = 'collision'
tallies[-1].scores = scores[do_nuclides]
if do_nuclides:
tallies[-1].nuclides = nuclides
tallies.append(openmc.Tally())
tallies[-1].filters = [mat_filter, energy_filter]
tallies[-1].estimator = 'tracklength'
tallies[-1].scores = scores[do_nuclides]
if do_nuclides:
tallies[-1].nuclides = nuclides
tallies.append(openmc.Tally())
tallies[-1].filters = [mat_filter, energy_filter,
energyout_filter]
tallies[-1].scores = ['scatter', 'nu-scatter', 'nu-fission']
if do_nuclides:
tallies[-1].nuclides = nuclides
self._input_set.tallies = openmc.Tallies(tallies)
self._input_set.export()
def _cleanup(self):
super(MGTalliesTestHarness, self)._cleanup()

2
tests/test_mgxs_library_ce_to_mg/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.140804E+00 2.937150E-02
1.140724E+00 2.952379E-02

7
tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
View file

@ -61,14 +61,14 @@ class MGXSTestHarness(PyAPITestHarness):
self._input_set.geometry = self.mgxs_lib.create_mg_mode()
# Modify settings so we can run in MG mode
self._input_set.settings.cross_sections = './mgxs.xml'
self._input_set.settings.cross_sections = './mgxs.h5'
self._input_set.settings.energy_mode = 'multi-group'
# Write modified input files
self._input_set.settings.export_to_xml()
self._input_set.geometry.export_to_xml()
self._input_set.materials.export_to_xml()
self._input_set.mgxs_file.export_to_xml()
self._input_set.mgxs_file.export_to_hdf5()
# Dont need tallies.xml, so remove the file
if os.path.exists('./tallies.xml'):
os.remove('./tallies.xml')
@ -83,9 +83,8 @@ class MGXSTestHarness(PyAPITestHarness):
returncode = openmc.run(openmc_exec=self._opts.exe)
def _cleanup(self):
return
super(MGXSTestHarness, self)._cleanup()
f = os.path.join(os.getcwd(), 'mgxs.xml')
f = os.path.join(os.getcwd(), 'mgxs.h5')
if os.path.exists(f):
os.remove(f)