From 976349cb516d0f950e9391bfe0f1039f30bcf9f4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Jul 2016 22:32:24 -0500 Subject: [PATCH] Calculate photon production yields properly for MF=13 data --- openmc/data/data.py | 18 ++++++++++++ openmc/data/function.py | 37 +++++++++++++++++++++++- openmc/data/neutron.py | 62 ++++++++++++++++++++++++++++++++++++++--- openmc/data/reaction.py | 39 +++++++++++++++----------- 4 files changed, 134 insertions(+), 22 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index 52043e3923..5f83d31c0c 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -156,3 +156,21 @@ REACTION_NAME.update({i: '(n,d{})'.format(i-650) for i in range(650, 699)}) REACTION_NAME.update({i: '(n,t{})'.format(i-700) for i in range(700, 749)}) REACTION_NAME.update({i: '(n,3He{})'.format(i-750) for i in range(750, 799)}) REACTION_NAME.update({i: '(n,a{})'.format(i-800) for i in range(800, 849)}) + +SUM_RULES = {1: [2, 3], + 3: [4, 5, 11, 16, 17, 22, 23, 24, 25, 27, 28, 29, 30, 32, 33, 34, 35, + 36, 37, 41, 42, 44, 45, 152, 153, 154, 156, 157, 158, 159, 160, + 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, + 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, + 186, 187, 188, 189, 190, 194, 195, 196, 198, 199, 200], + 4: list(range(50, 92)), + 16: list(range(875, 892)), + 18: [19, 20, 21, 38], + 27: [18, 101], + 101: [102, 103, 104, 105, 106, 107, 108, 109, 111, 112, 113, 114, + 115, 116, 117, 155, 182, 191, 192, 193, 197], + 103: list(range(600, 650)), + 104: list(range(650, 700)), + 105: list(range(700, 750)), + 106: list(range(750, 800)), + 107: list(range(800, 850))} diff --git a/openmc/data/function.py b/openmc/data/function.py index d9902794cf..e2efb4c101 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,4 +1,4 @@ -from collections import Iterable +from collections import Iterable, Callable from numbers import Real, Integral import numpy as np @@ -306,3 +306,38 @@ class Tabulated1D(object): y = ace.xss[idx + n_pairs:idx + 2*n_pairs] return Tabulated1D(x, y, breakpoints, interpolation) + + +class Sum(object): + """Sum of multiple functions. + + This class allows you to create a callable object which represents the sum + of other callable objects. This is used for summed reactions whereby the + cross section is defined as the sum of other cross sections. + + Parameters + ---------- + functions : Iterable of Callable + Functions which are to be added together + + Attributes + ---------- + functions : Iterable of Callable + Functions which are to be added together + + """ + + def __init__(self, functions): + self.functions = functions + + def __call__(self, x): + return sum(f(x) for f in self.functions) + + @property + def functions(self): + return self._functions + + @functions.setter + def functions(self, functions): + cv.check_type('functions', functions, Iterable, Callable) + self._functions = functions diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 30c1dbb5b1..c77ea86a1e 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -2,13 +2,14 @@ from __future__ import division, unicode_literals import sys from collections import OrderedDict, Iterable, Mapping from numbers import Integral, Real +from warnings import warn import numpy as np import h5py -from . import ATOMIC_SYMBOL +from .data import ATOMIC_SYMBOL, SUM_RULES from .ace import Table, get_table -from .function import Tabulated1D +from .function import Tabulated1D, Sum from .product import Product from .reaction import Reaction, _get_photon_products from .urr import ProbabilityTables @@ -84,6 +85,17 @@ class IncidentNeutron(object): self.summed_reactions = OrderedDict() self.urr = None + def __contains__(self, mt): + return mt in self.reactions or mt in self.summed_reactions + + def __getitem__(self, mt): + if mt in self.reactions: + return self.reactions[mt] + elif mt in self.summed_reactions: + return self.summed_reactions[mt] + else: + raise KeyError('No reaction with MT={}.'.format(mt)) + def __repr__(self): return "".format(self.name) @@ -190,6 +202,38 @@ class IncidentNeutron(object): (ProbabilityTables, type(None))) self._urr = urr + def get_reaction_components(self, mt): + """Determine what reactions make up summed reaction. + + Parameters + ---------- + mt : int + ENDF MT number of the reaction to find components of. + + Returns + ------- + mts : list of int + ENDF MT numbers of reactions that make up the summed reaction and + have cross sections provided. + + """ + if mt in self.reactions: + return [mt] + elif mt in SUM_RULES: + mts = SUM_RULES[mt] + complete = False + while not complete: + new_mts = [] + complete = True + for i, mt_i in enumerate(mts): + if mt_i in self.reactions: + new_mts.append(mt_i) + elif mt_i in SUM_RULES: + new_mts += SUM_RULES[mt_i] + complete = False + mts = new_mts + return mts + def export_to_hdf5(self, path, mode='a'): """Export table to an HDF5 file. @@ -386,9 +430,19 @@ class IncidentNeutron(object): n_photon_reactions].astype(int) for mt in np.unique(photon_mts // 1000): - if mt not in data.reactions: + if mt not in data: + if mt not in SUM_RULES: + warn('Photon production is present for MT={} but no ' + 'cross section is given.'.format(mt)) + continue + + # Create summed reaction with appropriate cross section rx = Reaction(mt) - rx.products += _get_photon_products(ace, mt) + mts = data.get_reaction_components(mt) + rx.xs = Sum([data.reactions[mt_i].xs for mt_i in mts]) + + # Determine summed cross section + rx.products += _get_photon_products(ace, rx) data.summed_reactions[mt] = rx # Read unresolved resonance probability tables diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 91d18bdce3..8fb79f109d 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1,5 +1,5 @@ from __future__ import division, unicode_literals -from collections import Iterable +from collections import Iterable, Callable from copy import deepcopy from numbers import Real from warnings import warn @@ -146,15 +146,15 @@ def _get_fission_products(ace): return products, derived_products -def _get_photon_products(ace, mt): +def _get_photon_products(ace, rx): """Generate photon products from an ACE table Parameters ---------- ace : openmc.data.ace.Table ACE table to read from - mt : int - MT number for the desired reaction + rx : openmc.data.Reaction + Reaction that generates photons Returns ------- @@ -176,9 +176,9 @@ def _get_photon_products(ace, mt): # MT=19,20,21,38, we assign the photon product to each of the individual # reactions if neutron_mt == 18: - if mt not in (18, 19, 20, 21, 38): + if rx.mt not in (18, 19, 20, 21, 38): continue - elif neutron_mt != mt: + elif neutron_mt != rx.mt: continue # Create photon product and assign to reactions @@ -205,15 +205,19 @@ def _get_photon_products(ace, mt): # Energy grid index at which data starts threshold_idx = int(ace.xss[idx]) - 1 - - # Get photon production cross section n_energy = int(ace.xss[idx + 1]) - photon._xs = ace.xss[idx + 2:idx + 2 + n_energy] - - # TODO: Determine yield based on ratio of cross sections energy = ace.xss[ace.jxs[1] + threshold_idx: ace.jxs[1] + threshold_idx + n_energy] - photon.yield_ = Tabulated1D(energy, photon._xs) + + # Get photon production cross section + photon_prod_xs = ace.xss[idx + 2:idx + 2 + n_energy] + neutron_xs = rx.xs(energy) + idx = np.where(neutron_xs > 0.) + + # Calculate photon yield + yield_ = np.zeros_like(photon_prod_xs) + yield_[idx] = photon_prod_xs[idx] / neutron_xs[idx] + photon.yield_ = Tabulated1D(energy, yield_) else: raise ValueError("MFTYPE must be 12, 13, 16. Got {0}".format( @@ -277,7 +281,7 @@ class Reaction(object): threshold_idx : int The index on the energy grid corresponding to the threshold of this reaction. - xs : openmc.data.Tabulated1D + xs : callable Microscopic cross section for this reaction as a function of incident energy products : Iterable of openmc.data.Product @@ -340,9 +344,10 @@ class Reaction(object): @xs.setter def xs(self, xs): - cv.check_type('reaction cross section', xs, Tabulated1D) - for y in xs.y: - cv.check_greater_than('reaction cross section', y, 0.0, True) + cv.check_type('reaction cross section', xs, Callable) + if isinstance(xs, Tabulated1D): + for y in xs.y: + cv.check_greater_than('reaction cross section', y, 0.0, True) self._xs = xs def to_hdf5(self, group): @@ -531,6 +536,6 @@ class Reaction(object): # ====================================================================== # PHOTON PRODUCTION - rx.products += _get_photon_products(ace, mt) + rx.products += _get_photon_products(ace, rx) return rx