mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Calculate photon production yields properly for MF=13 data
This commit is contained in:
Paul Romano
parent
20a2499587
commit
976349cb51
4 changed files with 134 additions and 22 deletions
|
|
@ -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))}
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 "<IncidentNeutron: {}>".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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue