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

Improve loading of cross sections from hdf5

Browse source 
Fix bug where reaction thresholds were not used.  Automatically build summed
reactions.  Allow Tab1 evaluations of un-ordered datasets.
This commit is contained in:
Sterling Harper 2016-07-27 18:02:38 -05:00
parent c8a188ca26
commit be7ca28dee
3 changed files with 34 additions and 27 deletions
Download patch file Download diff file Expand all files Collapse all files

48
openmc/data/function.py
View file

@ -80,50 +80,46 @@ class Tabulated1D(object):
# Get indices for interpolation
idx = np.searchsorted(self.x, x, side='right') - 1
# Find lowest valid index
i_low = np.searchsorted(idx, 0)
# Loop over interpolation regions
for k in range(len(self.breakpoints)):
# Determine which x values are within this interpolation range
i_high = np.searchsorted(idx, self.breakpoints[k] - 1)
# Get indices for the begining and ending of this region
i_begin = self.breakpoints[k-1] - 1 if k > 0 else 0
i_end = self.breakpoints[k] - 1
# Get x values and bounding (x,y) pairs
xk = x[i_low:i_high]
xi = self.x[idx[i_low:i_high]]
xi1 = self.x[idx[i_low:i_high] + 1]
yi = self.y[idx[i_low:i_high]]
yi1 = self.y[idx[i_low:i_high] + 1]
# Figure out which idx values lie within this region
contained = (idx >= i_begin) & (idx < i_end)
xk = x[contained] # x values in this region
xi = self.x[idx[contained]] # low edge of corresponding bins
xi1 = self.x[idx[contained] + 1] # high edge of corresponding bins
yi = self.y[idx[contained]]
yi1 = self.y[idx[contained] + 1]
if self.interpolation[k] == 1:
# Histogram
y[i_low:i_high] = yi
y[contined] = yi
elif self.interpolation[k] == 2:
# Linear-linear
y[i_low:i_high] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
y[contained] = yi + (xk - xi)/(xi1 - xi)*(yi1 - yi)
elif self.interpolation[k] == 3:
# Linear-log
y[i_low:i_high] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
y[contained] = yi + np.log(xk/xi)/np.log(xi1/xi)*(yi1 - yi)
elif self.interpolation[k] == 4:
# Log-linear
y[i_low:i_high] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
y[contained] = yi*np.exp((xk - xi)/(xi1 - xi)*np.log(yi1/yi))
elif self.interpolation[k] == 5:
# Log-log
y[i_low:i_high] = yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)*np.log(yi1/yi))
y[contained] = (yi*np.exp(np.log(xk/xi)/np.log(xi1/xi)
*np.log(yi1/yi)))
i_low = i_high
# In some cases, the first/last point of x may be less than the first
# value of self.x due only to precision, so we check if they're close
# and set them equal if so. Otherwise, the interpolated value might be
# out of range (and thus zero)
if np.isclose(x[0], self.x[0], 1e-8):
y[0] = self.y[0]
if np.isclose(x[-1], self.x[-1], 1e-8):
y[-1] = self.y[-1]
# In some cases, x values might be outside the tabulated region due only
# to precision, so we check if they're close and set them equal if so.
y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]
return y if iterable else y[0]

11
openmc/data/neutron.py
View file

@ -326,6 +326,17 @@ class IncidentNeutron(object):
tgroup = group['total_nu']
rx.derived_products.append(Product.from_hdf5(tgroup))
# Build summed reactions. Start from the highest MT number because high
# MTs never depend on lower MTs.
for mt_sum in sorted(SUM_RULES.keys())[::-1]:
if mt_sum not in data:
xs_components = [data[mt].xs for mt in SUM_RULES[mt_sum]
if mt in data]
if len(xs_components) > 0:
rxn = Reaction(mt_sum)
rxn.xs = Sum(xs_components)
data.summed_reactions[mt_sum] = rxn
# Read unresolved resonance probability tables
if 'urr' in group:
urr_group = group['urr']

2
openmc/data/reaction.py
View file

@ -400,7 +400,7 @@ class Reaction(object):
# Read cross section
if 'xs' in group:
xs = group['xs'].value
rx.xs = Tabulated1D(energy, xs)
rx.xs = Tabulated1D(energy[rx.threshold_idx:], xs)
# Determine number of products
n_product = 0