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

Rename openmc.data.container to openmc.data.function

Browse source
This commit is contained in:
Paul Romano 2016-07-19 21:41:41 -05:00
parent 5647db18b6
commit cf59208f34
10 changed files with 9 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

308
openmc/data/function.py Normal file
View file

@ -0,0 +1,308 @@
from collections import Iterable
from numbers import Real, Integral
import numpy as np
import openmc.checkvalue as cv
INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log',
4: 'log-linear', 5: 'log-log'}
class Tabulated1D(object):
"""A one-dimensional tabulated function.
This class mirrors the TAB1 type from the ENDF-6 format. A tabulated
function is specified by tabulated (x,y) pairs along with interpolation
rules that determine the values between tabulated pairs.
Once an object has been created, it can be used as though it were an actual
function, e.g.:
>>> f = Tabulated1D([0, 10], [4, 5])
>>> [f(xi) for xi in numpy.linspace(0, 10, 5)]
[4.0, 4.25, 4.5, 4.75, 5.0]
Parameters
----------
x : Iterable of float
Independent variable
y : Iterable of float
Dependent variable
breakpoints : Iterable of int
Breakpoints for interpolation regions
interpolation : Iterable of int
Interpolation scheme identification number, e.g., 3 means y is linear in
ln(x).
Attributes
----------
x : Iterable of float
Independent variable
y : Iterable of float
Dependent variable
breakpoints : Iterable of int
Breakpoints for interpolation regions
interpolation : Iterable of int
Interpolation scheme identification number, e.g., 3 means y is linear in
ln(x).
n_regions : int
Number of interpolation regions
n_pairs : int
Number of tabulated (x,y) pairs
"""
def __init__(self, x, y, breakpoints=None, interpolation=None):
if breakpoints is None or interpolation is None:
# Single linear-linear interpolation region by default
self.breakpoints = np.array([len(x)])
self.interpolation = np.array([2])
else:
self.breakpoints = np.asarray(breakpoints, dtype=int)
self.interpolation = np.asarray(interpolation, dtype=int)
self.x = np.asarray(x)
self.y = np.asarray(y)
def __call__(self, x):
# Check if input is array or scalar
if isinstance(x, Iterable):
iterable = True
x = np.array(x)
else:
iterable = False
x = np.array([x], dtype=float)
# Create output array
y = np.zeros_like(x)
# Get indices for interpolation
idx = np.searchsorted(self.x, x, side='right') - 1
# Find lowest valid index
i_low = np.searchsorted(idx, 0)
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 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]
if self.interpolation[k] == 1:
# Histogram
y[i_low:i_high] = yi
elif self.interpolation[k] == 2:
# Linear-linear
y[i_low:i_high] = 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)
elif self.interpolation[k] == 4:
# Log-linear
y[i_low:i_high] = 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))
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]
return y if iterable else y[0]
def __len__(self):
return len(self.x)
@property
def x(self):
return self._x
@property
def y(self):
return self._y
@property
def breakpoints(self):
return self._breakpoints
@property
def interpolation(self):
return self._interpolation
@property
def n_pairs(self):
return len(self.x)
@property
def n_regions(self):
return len(self.breakpoints)
@x.setter
def x(self, x):
cv.check_type('x values', x, Iterable, Real)
self._x = x
@y.setter
def y(self, y):
cv.check_type('y values', y, Iterable, Real)
self._y = y
@breakpoints.setter
def breakpoints(self, breakpoints):
cv.check_type('breakpoints', breakpoints, Iterable, Integral)
self._breakpoints = breakpoints
@interpolation.setter
def interpolation(self, interpolation):
cv.check_type('interpolation', interpolation, Iterable, Integral)
self._interpolation = interpolation
def integral(self):
"""Integral of the tabulated function over its tabulated range.
Returns
-------
numpy.ndarray
Array of same length as the tabulated data that represents partial
integrals from the bottom of the range to each tabulated point.
"""
# Create output array
partial_sum = np.zeros(len(self.x) - 1)
i_low = 0
for k in range(len(self.breakpoints)):
# Determine which x values are within this interpolation range
i_high = self.breakpoints[k] - 1
# Get x values and bounding (x,y) pairs
x0 = self.x[i_low:i_high]
x1 = self.x[i_low + 1:i_high + 1]
y0 = self.y[i_low:i_high]
y1 = self.y[i_low + 1:i_high + 1]
if self.interpolation[k] == 1:
# Histogram
partial_sum[i_low:i_high] = y0*(x1 - x0)
elif self.interpolation[k] == 2:
# Linear-linear
m = (y1 - y0)/(x1 - x0)
partial_sum[i_low:i_high] = (y0 - m*x0)*(x1 - x0) + \
m*(x1**2 - x0**2)/2
elif self.interpolation[k] == 3:
# Linear-log
logx = np.log(x1/x0)
m = (y1 - y0)/logx
partial_sum[i_low:i_high] = y0 + m*(x1*(logx - 1) + x0)
elif self.interpolation[k] == 4:
# Log-linear
m = np.log(y1/y0)/(x1 - x0)
partial_sum[i_low:i_high] = y0/m*(np.exp(m*(x1 - x0)) - 1)
elif self.interpolation[k] == 5:
# Log-log
m = np.log(y1/y0)/np.log(x1/x0)
partial_sum[i_low:i_high] = y0/((m + 1)*x0**m)*(
x1**(m + 1) - x0**(m + 1))
i_low = i_high
return np.concatenate(([0.], np.cumsum(partial_sum)))
def to_hdf5(self, group, name='xy'):
"""Write tabulated function to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : str
Name of the dataset to create
"""
dataset = group.create_dataset(name, data=np.vstack(
[self.x, self.y]))
dataset.attrs['type'] = np.string_('tab1')
dataset.attrs['breakpoints'] = self.breakpoints
dataset.attrs['interpolation'] = self.interpolation
@classmethod
def from_hdf5(cls, dataset):
"""Generate tabulated function from an HDF5 dataset
Parameters
----------
dataset : h5py.Dataset
Dataset to read from
Returns
-------
openmc.data.Tabulated1D
Function read from dataset
"""
x = dataset.value[0, :]
y = dataset.value[1, :]
breakpoints = dataset.attrs['breakpoints']
interpolation = dataset.attrs['interpolation']
return cls(x, y, breakpoints, interpolation)
@classmethod
def from_ace(cls, ace, idx=0):
"""Create a Tabulated1D object from an ACE table.
Parameters
----------
ace : openmc.data.ace.Table
An ACE table
idx : int
Offset to read from in XSS array (default of zero)
Returns
-------
openmc.data.Tabulated1D
Tabulated data object
"""
# Get number of regions and pairs
n_regions = int(ace.xss[idx])
n_pairs = int(ace.xss[idx + 1 + 2*n_regions])
# Get interpolation information
idx += 1
if n_regions > 0:
breakpoints = ace.xss[idx:idx + n_regions].astype(int)
interpolation = ace.xss[idx + n_regions:idx + 2*n_regions].astype(int)
else:
# 0 regions implies linear-linear interpolation by default
breakpoints = np.array([n_pairs])
interpolation = np.array([2])
# Get (x,y) pairs
idx += 2*n_regions + 1
x = ace.xss[idx:idx + n_pairs]
y = ace.xss[idx + n_pairs:idx + 2*n_pairs]
return Tabulated1D(x, y, breakpoints, interpolation)