diff --git a/docs/source/_static/theme_overrides.css b/docs/source/_static/theme_overrides.css new file mode 100644 index 000000000..7c1a52022 --- /dev/null +++ b/docs/source/_static/theme_overrides.css @@ -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; +} diff --git a/docs/source/conf.py b/docs/source/conf.py index 32154fa60..65db07b25 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -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. diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 2553e856b..f6ac53482 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1475,114 +1475,203 @@ The ```` 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| + | |``scatter-2``. | + +----------------------+---------------------------------------------------+ + |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 `` scatter-P2 ``. | + +----------------------+---------------------------------------------------+ + |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 `` scatter-2 - ``. 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 `` - scatter-P2 ``. 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. diff --git a/openmc/aggregate.py b/openmc/aggregate.py deleted file mode 100644 index 011ad3850..000000000 --- a/openmc/aggregate.py +++ /dev/null @@ -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 diff --git a/openmc/cross.py b/openmc/arithmetic.py similarity index 52% rename from openmc/cross.py rename to openmc/arithmetic.py index ee0fcb4c8..65e12c5d2 100644 --- a/openmc/cross.py +++ b/openmc/arithmetic.py @@ -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 diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index ca9fa5d66..4b8d8a914 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -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) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 0eec6d458..f24127458 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -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') diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 42e957506..ab37be78d 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -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) diff --git a/openmc/tallies.py b/openmc/tallies.py index d202497aa..da209f527 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -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: diff --git a/src/ace.F90 b/src/ace.F90 index e97338b55..8d73f8bd4 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -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. diff --git a/src/ace_header.F90 b/src/ace_header.F90 index d7c298979..2331f9117 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -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) diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 new file mode 100644 index 000000000..282d51b73 --- /dev/null +++ b/src/angle_distribution.F90 @@ -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 diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 index b3d77e2f8..f596536d0 100644 --- a/src/distribution_univariate.F90 +++ b/src/distribution_univariate.F90 @@ -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 diff --git a/src/endf.F90 b/src/endf.F90 index ba324722c..64f26539a 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -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) diff --git a/src/endf_header.F90 b/src/endf_header.F90 index af62231a5..7388ea2f5 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -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 diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 new file mode 100644 index 000000000..db3dcc441 --- /dev/null +++ b/src/energy_distribution.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index ebbdb58b2..8eb70d2e6 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -23,7 +23,6 @@ module global #endif implicit none - save ! ============================================================================ ! GEOMETRY-RELATED VARIABLES diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f811cab04..0a90b5317 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 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 diff --git a/src/interpolation.F90 b/src/interpolation.F90 index 9e28bc086..5f8787067 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -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 diff --git a/src/output.F90 b/src/output.F90 index e09797b28..f9b025f4e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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_,*) diff --git a/src/physics.F90 b/src/physics.F90 index 217d23218..31eb981ab 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,6 +1,6 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy + use ace_header, only: Nuclide, Reaction use constants use cross_section, only: elastic_xs_0K use endf, only: reaction_name @@ -17,12 +17,10 @@ module physics use random_lcg, only: prn use search, only: binary_search use string, only: to_str + use secondary_uncorrelated, only: UncorrelatedAngleEnergy implicit none -! TODO: Figure out how to write particle restart files in sample_angle, -! sample_energy, etc. - contains !=============================================================================== @@ -458,7 +456,10 @@ contains vel = sqrt(dot_product(v_n, v_n)) ! Sample scattering angle - mu_cm = sample_angle(rxn, E) + select type (dist => rxn%secondary%distribution(1)%obj) + type is (UncorrelatedAngleEnergy) + mu_cm = dist%angle%sample(E) + end select ! Determine direction cosines in CM uvw_cm = v_n/vel @@ -1195,16 +1196,13 @@ contains integer :: lc ! index before start of energies/nu values integer :: NR ! number of interpolation regions integer :: NE ! number of energies tabulated - integer :: law ! energy distribution law integer :: n_sample ! number of times resampling real(8) :: nu_t ! total nu real(8) :: nu_d ! delayed nu - real(8) :: mu ! fission neutron angular cosine real(8) :: beta ! delayed neutron fraction real(8) :: xi ! random number real(8) :: yield ! delayed neutron precursor yield real(8) :: prob ! cumulative probability - type(DistEnergy), pointer :: edist ! Determine total nu nu_t = nu_total(nuc, p % E) @@ -1248,18 +1246,13 @@ contains ! set the delayed group for the particle born from fission p % delayed_group = j - ! select energy distribution for group j - law = nuc % nu_d_edist(j) % law - edist => nuc % nu_d_edist(j) - ! sample from energy distribution n_sample = 0 do - if (law == 44 .or. law == 61) then - call sample_energy(edist, p % E, E_out, mu) - else - call sample_energy(edist, p % E, E_out) - end if + select type (aedist => nuc%nu_d_edist(j)%obj) + type is (UncorrelatedAngleEnergy) + E_out = aedist%energy%sample(p%E) + end select ! resample if energy is greater than maximum neutron energy if (E_out < energy_max_neutron) exit @@ -1281,14 +1274,9 @@ contains p % delayed_group = 0 ! sample from prompt neutron energy distribution - law = rxn % edist % law n_sample = 0 do - if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, p % E, E_out, prob) - else - call sample_energy(rxn%edist, p % E, E_out) - end if + call rxn%secondary%sample(p%E, E_out, prob) ! resample if energy is greater than maximum neutron energy if (E_out < energy_max_neutron) exit @@ -1317,41 +1305,18 @@ contains type(Particle), intent(inout) :: p integer :: i ! loop index - integer :: law ! secondary energy distribution law real(8) :: E ! energy in lab (incoming/outgoing) real(8) :: mu ! cosine of scattering angle in lab real(8) :: A ! atomic weight ratio of nuclide real(8) :: E_in ! incoming energy real(8) :: E_cm ! outgoing energy in center-of-mass - real(8) :: Q ! Q-value of reaction real(8) :: yield ! neutron yield ! copy energy of neutron E_in = p % E - ! determine A and Q - A = nuc % awr - Q = rxn % Q_value - - ! determine secondary energy distribution law - law = rxn % edist % law - - ! sample scattering angle - mu = sample_angle(rxn, E_in) - - ! sample outgoing energy - if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, E_in, E, mu) - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) - elseif (law == 66) then - call sample_energy(rxn%edist, E_in, E, A=A, Q=Q) - else - call sample_energy(rxn%edist, E_in, E) - end if + ! sample outgoing energy and scattering cosine + call rxn%secondary%sample(E_in, E, mu) ! if scattering system is in center-of-mass, transfer cosine of scattering ! angle and outgoing energy from CM to LAB @@ -1359,19 +1324,19 @@ contains E_cm = E ! determine outgoing energy in lab + A = nuc%awr E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) end if + ! Because of floating-point roundoff, it may be possible for mu to be + ! outside of the range [-1,1). In these cases, we just set mu to exactly -1 + ! or 1 + if (abs(mu) > ONE) mu = sign(ONE,mu) + ! Set outgoing energy and scattering angle p % E = E p % mu = mu @@ -1391,893 +1356,4 @@ contains end subroutine inelastic_scatter -!=============================================================================== -! SAMPLE_ANGLE samples the cosine of the angle between incident and exiting -! particle directions either from 32 equiprobable bins or from a tabular -! distribution. -!=============================================================================== - - function sample_angle(rxn, E) result(mu) - type(Reaction), intent(in) :: rxn ! reaction - real(8), intent(in) :: E ! incoming energy - - real(8) :: xi ! random number on [0,1) - integer :: interp ! type of interpolation - integer :: type ! angular distribution type - integer :: i ! incoming energy bin - integer :: n ! number of incoming energy bins - integer :: lc ! location in data array - integer :: NP ! number of points in cos distribution - integer :: k ! index on cosine grid - real(8) :: r ! interpolation factor on incoming energy - real(8) :: frac ! interpolation fraction on cosine - real(8) :: mu0 ! cosine in bin k - real(8) :: mu1 ! cosine in bin k+1 - real(8) :: mu ! final cosine sampled - real(8) :: c_k ! cumulative frequency at k - real(8) :: c_k1 ! cumulative frequency at k+1 - real(8) :: p0,p1 ! probability distribution - - ! check if reaction has angular distribution -- if not, sample outgoing - ! angle isotropically - if (.not. rxn % has_angle_dist) then - mu = TWO * prn() - ONE - return - end if - - ! determine number of incoming energies - n = rxn % adist % n_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 < rxn % adist % energy(1)) then - i = 1 - r = ZERO - elseif (E > rxn % adist % energy(n)) then - i = n - 1 - r = ONE - else - i = binary_search(rxn % adist % energy, n, E) - r = (E - rxn % adist % energy(i)) / & - (rxn % adist % energy(i+1) - rxn % adist % energy(i)) - end if - - ! Sample between the ith and (i+1)th bin - if (r > prn()) i = i + 1 - - ! check whether this is a 32-equiprobable bin or a tabular distribution - lc = rxn % adist % location(i) - type = rxn % adist % type(i) - if (type == ANGLE_ISOTROPIC) then - mu = TWO * prn() - ONE - elseif (type == ANGLE_32_EQUI) then - ! sample cosine bin - xi = prn() - k = 1 + int(32.0_8*xi) - - ! calculate cosine - mu0 = rxn % adist % data(lc + k) - mu1 = rxn % adist % data(lc + k+1) - mu = mu0 + (32.0_8 * xi - k + ONE) * (mu1 - mu0) - - elseif (type == ANGLE_TABULAR) then - interp = int(rxn % adist % data(lc + 1)) - NP = int(rxn % adist % data(lc + 2)) - - ! determine outgoing cosine bin - xi = prn() - lc = lc + 2 - c_k = rxn % adist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = rxn % adist % data(lc + 2*NP + k+1) - if (xi < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - p0 = rxn % adist % data(lc + NP + k) - mu0 = rxn % adist % data(lc + k) - if (interp == HISTOGRAM) then - ! Histogram interpolation - if (p0 > ZERO) then - mu = mu0 + (xi - c_k)/p0 - else - mu = mu0 - end if - - elseif (interp == LINEAR_LINEAR) then - ! Linear-linear interpolation - p1 = rxn % adist % data(lc + NP + k+1) - mu1 = rxn % adist % data(lc + k+1) - - frac = (p1 - p0)/(mu1 - mu0) - if (frac == ZERO) then - mu = mu0 + (xi - c_k)/p0 - else - mu = mu0 + (sqrt(max(ZERO, p0*p0 + 2*frac*(xi - c_k))) - p0)/frac - end if - else - ! call write_particle_restart(p) - call fatal_error("Unknown interpolation type: " // trim(to_str(interp))) - end if - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) - - else - ! call write_particle_restart(p) - call fatal_error("Unknown angular distribution type: " & - &// trim(to_str(type))) - end if - - end function sample_angle - -!=============================================================================== -! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary -! neutron from a collision or for a prompt/delayed fission neutron -!=============================================================================== - - recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q) - type(DistEnergy), intent(in) :: edist - real(8), intent(in) :: E_in ! incoming energy of neutron - real(8), intent(out) :: E_out ! outgoing energy - real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle - real(8), intent(in), optional :: A ! mass number of nuclide - real(8), intent(in), optional :: Q ! Q-value of reaction - - integer :: i ! index on incoming energy grid - integer :: k ! sampled index on outgoing grid - integer :: l ! sampled index on incoming grid - integer :: n_sample ! number of rejections - integer :: lc ! dummy index - integer :: NR ! number of interpolation regions - integer :: NE ! number of energies - integer :: NET ! number of outgoing energies - integer :: INTTp ! combination of INTT and ND - integer :: INTT ! 1 = histogram, 2 = linear-linear - integer :: JJ ! 1 = histogram, 2 = linear-linear - integer :: ND ! number of discrete lines - integer :: NP ! number of points in distribution - - real(8) :: p_valid ! probability of law validity - - 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_A ! Kalbach-Mann parameter R - real(8) :: KM_R ! Kalbach-Mann parameter R - real(8) :: A_k, A_k1 ! Kalbach-Mann A on outgoing grid l - real(8) :: R_k, R_k1 ! Kalbach-Mann R on outgoing grid l - - real(8) :: Watt_a, Watt_b ! Watt spectrum parameters - - real(8) :: mu_k ! angular cosine in bin k - real(8) :: mu_k1 ! angular cosine in bin k+1 - real(8) :: p_k ! angular pdf in bin k - real(8) :: p_k1 ! angular pdf in bin k+1 - - real(8) :: r ! interpolation factor on incoming energy - real(8) :: frac ! interpolation factor on outgoing energy - real(8) :: U ! restriction energy - real(8) :: T ! nuclear temperature - - real(8) :: Ap ! total mass ratio for n-body dist - integer :: n_bodies ! number of bodies for n-body dist - real(8) :: E_max ! parameter for n-body dist - real(8) :: x, y, v ! intermediate variables for n-body dist - real(8) :: r1, r2, r3, r4, r5, r6 - logical :: histogram_interp ! use histogram interpolation on incoming energy - - ! ========================================================================== - ! SAMPLE ENERGY DISTRIBUTION IF THERE ARE MULTIPLE - - if (associated(edist % next)) then - p_valid = interpolate_tab1(edist % p_valid, E_in) - - if (prn() > p_valid) then - if (edist % law == 44 .or. edist % law == 61) then - call sample_energy(edist%next, E_in, E_out, mu_out) - elseif (edist % law == 66) then - call sample_energy(edist%next, E_in, E_out, A=A, Q=Q) - else - call sample_energy(edist%next, E_in, E_out) - end if - return - end if - end if - - ! Determine which secondary energy distribution law to use - select case (edist % law) - case (1) - ! ======================================================================= - ! TABULAR EQUIPROBABLE ENERGY BINS - - ! read number of interpolation regions, incoming energies, and outgoing - ! energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - NET = int(edist % data(3 + 2*NR + NE)) - if (NR > 0) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample equiprobable energy bins.") - end if - - ! determine index on incoming energy grid and interpolation factor - lc = 2 + 2*NR - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - - ! Sample outgoing energy bin - r1 = prn() - k = 1 + int(NET * r1) - - ! Determine E_1 and E_K - lc = 3 + 3*NR + NE + (i-1)*NET - E_i_1 = edist % data(lc + 1) - E_i_K = edist % data(lc + NET) - - lc = 3 + 3*NR + NE + i*NET - E_i1_1 = edist % data(lc + 1) - E_i1_K = edist % data(lc + NET) - - 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 - lc = 3 + 2*NR + NE + (l-1)*NET - E_l_k = edist % data(lc+k) - E_l_k1 = edist % data(lc+k+1) - - ! Determine E' (denoted here as E_out) - r2 = prn() - E_out = E_l_k + r2*(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 - - case (3) - ! ======================================================================= - ! INELASTIC LEVEL SCATTERING - - E_out = edist%data(2) * (E_in - edist%data(1)) - - case (4) - ! ======================================================================= - ! CONTINUOUS TABULAR DISTRIBUTION - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR == 1) then - histogram_interp = (edist % data(3) == 1) - else if (NR > 1) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample continuous tabular distribution.") - 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 - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - if (histogram_interp) then - l = i - else - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - end if - - ! interpolation for energy E1 and EK - lc = int(edist%data(2 + 2*NR + NE + i)) - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - lc = int(edist%data(2 + 2*NR + NE + i + 1)) - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - 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 location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart(p) - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == 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 (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+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 + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - else - ! call write_particle_restart(p) - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - 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 - - case (5) - ! ======================================================================= - ! GENERAL EVAPORATION SPECTRUM - - case (7) - ! ======================================================================= - ! MAXWELL FISSION SPECTRUM - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine nuclear temperature from tabulated function - T = interpolate_tab1(edist % data, E_in) - - ! determine restriction energy - lc = 2 + 2*NR + 2*NE - U = edist % data(lc + 1) - - n_sample = 0 - do - ! sample maxwell fission spectrum - E_out = maxwell_spectrum(T) - - ! accept energy based on restriction energy - if (E_out <= E_in - U) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on Maxwell fission spectrum.") - end if - end do - - case (9) - ! ======================================================================= - ! EVAPORATION SPECTRUM - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine nuclear temperature from tabulated function - T = interpolate_tab1(edist % data, E_in) - - ! determine restriction energy - lc = 2 + 2*NR + 2*NE - U = edist % data(lc + 1) - - y = (E_in - U)/T - v = 1 - exp(-y) - - ! sample outgoing energy based on evaporation spectrum probability - ! density function - n_sample = 0 - do - x = -log((1 - v*prn())*(1 - v*prn())) - if (x <= y) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on evaporation spectrum.") - end if - end do - - E_out = x*T - - case (11) - ! ======================================================================= - ! ENERGY-DEPENDENT WATT SPECTRUM - - ! read number of interpolation regions and incoming energies for - ! parameter 'a' - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - - ! determine Watt parameter 'a' from tabulated function - Watt_a = interpolate_tab1(edist % data, E_in) - - ! determine Watt parameter 'b' from tabulated function - lc = 2 + 2*(NR + NE) - Watt_b = interpolate_tab1(edist % data, E_in, lc + 1) - - ! read number of interpolation regions and incoming energies for - ! parameter 'a' - NR = int(edist % data(lc + 1)) - NE = int(edist % data(lc + 2 + 2*NR)) - - ! determine restriction energy - lc = lc + 2 + 2*(NR + NE) - U = edist % data(lc + 1) - - n_sample = 0 - do - ! Sample energy-dependent Watt fission spectrum - E_out = watt_spectrum(Watt_a, Watt_b) - - ! accept energy based on restriction energy - if (E_out <= E_in - U) exit - - ! check for large number of rejections - n_sample = n_sample + 1 - if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) - call fatal_error("Too many rejections on Watt spectrum.") - end if - end do - - case (44) - ! ======================================================================= - ! KALBACH-MANN CORRELATED SCATTERING - - if (.not. present(mu_out)) then - ! call write_particle_restart(p) - call fatal_error("Law 44 called without giving mu_out as argument.") - end if - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR > 0) then - ! call write_particle_restart(p) - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample Kalbach-Mann distribution.") - 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 - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - - ! determine endpoints on grid i - lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - ! determine endpoints on grid i+1 - lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1 - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - 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 location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart(p) - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == 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 = edist % data(lc + 3*NP + k) - KM_A = edist % data(lc + 4*NP + k) - - elseif (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+k+1) - - ! Find E prime - 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 + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - - ! Determine Kalbach-Mann parameters - R_k = edist % data(lc + 3*NP + k) - R_k1 = edist % data(lc + 3*NP + k+1) - A_k = edist % data(lc + 4*NP + k) - A_k1 = edist % data(lc + 4*NP + k+1) - - KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) - KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - 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 - r3 = prn() - r4 = prn() - if (r3 > KM_R) then - T = (TWO*r4 - ONE) * sinh(KM_A) - mu_out = log(T + sqrt(T*T + ONE))/KM_A - else - mu_out = log(r4*exp(KM_A) + (ONE - r4)*exp(-KM_A))/KM_A - end if - - case (61) - ! ======================================================================= - ! CORRELATED ENERGY AND ANGLE DISTRIBUTION - - if (.not. present(mu_out)) then - ! call write_particle_restart() - call fatal_error("Law 61 called without giving mu_out as argument.") - end if - - ! read number of interpolation regions and incoming energies - NR = int(edist % data(1)) - NE = int(edist % data(2 + 2*NR)) - if (NR > 0) then - ! call write_particle_restart() - call fatal_error("Multiple interpolation regions not supported while & - &attempting to sample correlated energy-angle distribution.") - 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 - lc = 2 + 2*NR - if (E_in < edist % data(lc+1)) then - i = 1 - r = ZERO - elseif (E_in > edist % data(lc+NE)) then - i = NE - 1 - r = ONE - else - i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & - (edist%data(lc+i+1) - edist%data(lc+i)) - end if - - ! Sample between the ith and (i+1)th bin - r2 = prn() - if (r > r2) then - l = i + 1 - else - l = i - end if - - ! determine endpoints on grid i - lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i - NP = int(edist%data(lc + 2)) - E_i_1 = edist%data(lc + 2 + 1) - E_i_K = edist%data(lc + 2 + NP) - - ! determine endpoints on grid i+1 - lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1 - NP = int(edist%data(lc + 2)) - E_i1_1 = edist%data(lc + 2 + 1) - E_i1_K = edist%data(lc + 2 + NP) - - 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 location of outgoing energies, pdf, cdf for E(l) - lc = int(edist % data(2 + 2*NR + NE + l)) - - ! determine type of interpolation and number of discrete lines - INTTp = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - if (INTTp > 10) then - INTT = mod(INTTp,10) - ND = (INTTp - INTT)/10 - else - INTT = INTTp - ND = 0 - end if - - if (ND > 0) then - ! discrete lines present - ! call write_particle_restart() - call fatal_error("Discrete lines in continuous tabular distributed not & - &yet supported") - end if - - ! determine outgoing energy bin - r1 = prn() - lc = lc + 2 ! start of EOUT - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - E_l_k = edist % data(lc+k) - p_l_k = edist % data(lc+NP+k) - if (INTT == 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 (INTT == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = edist % data(lc+k+1) - p_l_k1 = edist % data(lc+NP+k+1) - - ! Find E prime - 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 + & - 2*frac*(r1 - c_k))) - p_l_k)/frac - end if - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(INTT))) - 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 - lc = int(edist % data(lc + 3*NP + k)) - else - lc = int(edist % data(lc + 3*NP + k + 1)) - end if - - ! Check if angular distribution is isotropic - if (lc == 0) then - mu_out = TWO * prn() - ONE - return - end if - - ! interpolation type and number of points in angular distribution - JJ = int(edist % data(lc + 1)) - NP = int(edist % data(lc + 2)) - - ! determine outgoing cosine bin - r3 = prn() - lc = lc + 2 - c_k = edist % data(lc + 2*NP + 1) - do k = 1, NP - 1 - c_k1 = edist % data(lc + 2*NP + k+1) - if (r3 < c_k1) exit - c_k = c_k1 - end do - - ! check to make sure k is <= NP - 1 - k = min(k, NP - 1) - - p_k = edist % data(lc + NP + k) - mu_k = edist % data(lc + k) - if (JJ == HISTOGRAM) then - ! Histogram interpolation - if (p_k > ZERO) then - mu_out = mu_k + (r3 - c_k)/p_k - else - mu_out = mu_k - end if - - elseif (JJ == LINEAR_LINEAR) then - ! Linear-linear interpolation - p_k1 = edist % data(lc + NP + k+1) - mu_k1 = edist % data(lc + k+1) - - frac = (p_k1 - p_k)/(mu_k1 - mu_k) - if (frac == ZERO) then - mu_out = mu_k + (r3 - c_k)/p_k - else - mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac - end if - else - ! call write_particle_restart() - call fatal_error("Unknown interpolation type: " // trim(to_str(JJ))) - end if - - case (66) - ! ======================================================================= - ! N-BODY PHASE SPACE DISTRIBUTION - - ! read number of bodies in phase space and total mass ratio - n_bodies = int(edist % data(1)) - Ap = edist % data(2) - - ! determine E_max parameter - E_max = (Ap - ONE)/Ap * (A/(A+ONE)*E_in + Q) - - ! x is essentially a Maxwellian distribution - x = maxwell_spectrum(ONE) - - select case (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 - - case (67) - ! ======================================================================= - ! LABORATORY ENERGY-ANGLE LAW - - end select - - end subroutine sample_energy - end module physics diff --git a/src/search.F90 b/src/search.F90 index 0c345471c..f338105f4 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -1,7 +1,6 @@ module search use constants - use error, only: fatal_error implicit none diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90 new file mode 100644 index 000000000..c0289d55e --- /dev/null +++ b/src/secondary_correlated.F90 @@ -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 diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 new file mode 100644 index 000000000..d85807412 --- /dev/null +++ b/src/secondary_header.F90 @@ -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 diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 new file mode 100644 index 000000000..5e6949206 --- /dev/null +++ b/src/secondary_kalbach.F90 @@ -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 diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 new file mode 100644 index 000000000..22a56aa12 --- /dev/null +++ b/src/secondary_uncorrelated.F90 @@ -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 diff --git a/src/state_point.F90 b/src/state_point.F90 index 4c9d6f000..7991719c8 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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) diff --git a/src/summary.F90 b/src/summary.F90 index 85d014692..c31068eb9 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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) diff --git a/tests/test_energy_laws/geometry.xml b/tests/test_energy_laws/geometry.xml new file mode 100644 index 000000000..c42f45597 --- /dev/null +++ b/tests/test_energy_laws/geometry.xml @@ -0,0 +1,5 @@ + + + + + diff --git a/tests/test_energy_laws/materials.xml b/tests/test_energy_laws/materials.xml new file mode 100644 index 000000000..97f9b84be --- /dev/null +++ b/tests/test_energy_laws/materials.xml @@ -0,0 +1,11 @@ + + + 71c + + + + + + + + diff --git a/tests/test_energy_laws/results_true.dat b/tests/test_energy_laws/results_true.dat new file mode 100644 index 000000000..48eb6bc81 --- /dev/null +++ b/tests/test_energy_laws/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +2.130076E+00 1.938907E-03 diff --git a/tests/test_energy_laws/settings.xml b/tests/test_energy_laws/settings.xml new file mode 100644 index 000000000..1c3f444f0 --- /dev/null +++ b/tests/test_energy_laws/settings.xml @@ -0,0 +1,11 @@ + + + + 10 + 5 + 1000 + + + + + diff --git a/tests/test_energy_laws/test_energy_laws.py b/tests/test_energy_laws/test_energy_laws.py new file mode 100644 index 000000000..0e593876f --- /dev/null +++ b/tests/test_energy_laws/test_energy_laws.py @@ -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()