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

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

Browse source
This commit is contained in:
Sterling Harper 2016-01-30 23:19:16 -05:00
commit fe6f24307a
33 changed files with 2226 additions and 2042 deletions
Download patch file Download diff file Expand all files Collapse all files

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

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

4
docs/source/conf.py
View file

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

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

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

409
openmc/aggregate.py
View file

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

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

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

2
openmc/mgxs/library.py
View file

@ -396,7 +396,7 @@ class Library(object):
cv.check_type('statepoint', statepoint, openmc.StatePoint)
if not statepoint.with_summary:
if statepoint.summary is None:
msg = 'Unable to load data from a statepoint which has not been ' \
'linked with a summary file'
raise ValueError(msg)

39
openmc/mgxs/mgxs.py
View file

@ -583,7 +583,7 @@ class MGXS(object):
cv.check_type('statepoint', statepoint, openmc.statepoint.StatePoint)
if not statepoint.with_summary:
if statepoint.summary is None:
msg = 'Unable to load data from a statepoint which has not been ' \
'linked with a summary file'
raise ValueError(msg)
@ -612,6 +612,11 @@ class MGXS(object):
filters = []
filter_bins = []
# Clear any tallies previously loaded from a statepoint
self._tallies = None
self._xs_tally = None
self._rxn_rate_tally = None
# Find, slice and store Tallies from StatePoint
# The tally slicing is needed if tally merging was used
for tally_type, tally in self.tallies.items():
@ -673,14 +678,14 @@ class MGXS(object):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append('energy')
@ -838,7 +843,7 @@ class MGXS(object):
"""
# Construct a collection of the subdomain filter bins to average across
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains)
@ -881,7 +886,7 @@ class MGXS(object):
"""
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1010,7 +1015,7 @@ class MGXS(object):
xs_results = h5py.File(filename, 'w')
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1192,7 +1197,7 @@ class MGXS(object):
"""
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
if nuclides != 'all' and nuclides != 'sum':
cv.check_iterable_type('nuclides', nuclides, basestring)
@ -1252,7 +1257,7 @@ class MGXS(object):
columns = ['group in']
# Select out those groups the user requested
if groups != 'all':
if not isinstance(groups, basestring):
if 'group in' in df:
df = df[df['group in'].isin(groups)]
if 'group out' in df:
@ -1789,21 +1794,21 @@ class ScatterMatrixXS(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if in_groups != 'all':
if not isinstance(in_groups, basestring):
cv.check_iterable_type('groups', in_groups, Integral)
for group in in_groups:
filters.append('energy')
filter_bins.append((self.energy_groups.get_group_bounds(group),))
# Construct list of energy group bounds tuples for all requested groups
if out_groups != 'all':
if not isinstance(out_groups, basestring):
cv.check_iterable_type('groups', out_groups, Integral)
for group in out_groups:
filters.append('energyout')
@ -1887,7 +1892,7 @@ class ScatterMatrixXS(MGXS):
"""
# Construct a collection of the subdomains to report
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
elif self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
@ -1926,12 +1931,6 @@ class ScatterMatrixXS(MGXS):
bounds = self.energy_groups.get_group_bounds(group)
string += template.format('', group, bounds[0], bounds[1])
if subdomains == 'all':
if self.domain_type == 'distribcell':
subdomains = np.arange(self.num_subdomains, dtype=np.int)
else:
subdomains = [self.domain.id]
# Loop over all subdomains
for subdomain in subdomains:
@ -2135,14 +2134,14 @@ class Chi(MGXS):
filter_bins = []
# Construct a collection of the domain filter bins
if subdomains != 'all':
if not isinstance(subdomains, basestring):
cv.check_iterable_type('subdomains', subdomains, Integral)
for subdomain in subdomains:
filters.append(self.domain_type)
filter_bins.append((subdomain,))
# Construct list of energy group bounds tuples for all requested groups
if groups != 'all':
if not isinstance(groups, basestring):
cv.check_iterable_type('groups', groups, Integral)
for group in groups:
filters.append('energyout')

4
openmc/statepoint.py
View file

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

11
openmc/tallies.py
View file

@ -12,8 +12,7 @@ import sys
import numpy as np
from openmc import Mesh, Filter, Trigger, Nuclide, TallyDerivative
from openmc.cross import CrossScore, CrossNuclide, CrossFilter
from openmc.aggregate import AggregateScore, AggregateNuclide, AggregateFilter
from openmc.arithmetic import *
from openmc.filter import _FILTER_TYPES
import openmc.checkvalue as cv
from openmc.clean_xml import *
@ -2717,11 +2716,11 @@ class Tally(object):
new_tally = copy.deepcopy(self)
new_tally.sparse = False
if self.sum is not None:
if not self.derived and self.sum is not None:
new_sum = self.get_values(scores, filters, filter_bins,
nuclides, 'sum')
new_tally.sum = new_sum
if self.sum_sq is not None:
if not self.derived and self.sum_sq is not None:
new_sum_sq = self.get_values(scores, filters, filter_bins,
nuclides, 'sum_sq')
new_tally.sum_sq = new_sum_sq
@ -2992,10 +2991,10 @@ class Tally(object):
diag_indices[start:end] = indices + (i * new_filter.num_bins**2)
# Inject this Tally's data along the diagonal of the diagonalized Tally
if self.sum is not None:
if not self.derived and self.sum is not None:
new_tally._sum = np.zeros(new_tally.shape, dtype=np.float64)
new_tally._sum[diag_indices, :, :] = self.sum
if self.sum_sq is not None:
if not self.derived and self.sum_sq is not None:
new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64)
new_tally._sum_sq[diag_indices, :, :] = self.sum_sq
if self.mean is not None:

770
src/ace.F90
View file

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

77
src/ace_header.F90
View file

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

63
src/angle_distribution.F90 Normal file
View file

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

29
src/distribution_univariate.F90
View file

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

47
src/endf.F90
View file

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

34
src/endf_header.F90
View file

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

429
src/energy_distribution.F90 Normal file
View file

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

1
src/global.F90
View file

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

31
src/input_xml.F90
View file

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

1
src/interpolation.F90
View file

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

40
src/output.F90
View file

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

962
src/physics.F90
View file

File diff suppressed because it is too large Load diff

1
src/search.F90
View file

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

148
src/secondary_correlated.F90 Normal file
View file

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

78
src/secondary_header.F90 Normal file
View file

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

164
src/secondary_kalbach.F90 Normal file
View file

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

48
src/secondary_uncorrelated.F90 Normal file
View file

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

49
src/state_point.F90
View file

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

49
src/summary.F90
View file

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

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

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

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

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

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

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

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

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

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

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