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

Extend openmc.data to be able to import ENDF data and perform resonance

Browse source 
reconstruction.
This commit is contained in:
Paul Romano 2016-02-24 11:01:12 -06:00
parent 4f95b4a7f6
commit cbcd3a7f57
21 changed files with 3308 additions and 171 deletions
Download patch file Download diff file Expand all files Collapse all files

5
.gitignore vendored
View file

@ -89,3 +89,8 @@ docs/source/pythonapi/examples/mgxs
docs/source/pythonapi/examples/tracks
docs/source/pythonapi/examples/fission-rates
docs/source/pythonapi/examples/plots
# Cython files
*.c
*.html
*.so

73
docs/source/pythonapi/index.rst
View file

@ -333,6 +333,7 @@ Functions
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.model.create_triso_lattice
openmc.model.pack_trisos
@ -341,16 +342,6 @@ Functions
:mod:`openmc.data` -- Nuclear Data Interface
--------------------------------------------
Physical Data
-------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.data.atomic_mass
Core Classes
------------
@ -363,9 +354,20 @@ Core Classes
openmc.data.Reaction
openmc.data.Product
openmc.data.Tabulated1D
openmc.data.FissionEnergyRelease
openmc.data.ThermalScattering
openmc.data.CoherentElastic
openmc.data.FissionEnergyRelease
Core Functions
--------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.data.atomic_mass
openmc.data.write_compact_458_library
Angle-Energy Distributions
--------------------------
@ -380,6 +382,7 @@ Angle-Energy Distributions
openmc.data.CorrelatedAngleEnergy
openmc.data.UncorrelatedAngleEnergy
openmc.data.NBodyPhaseSpace
openmc.data.LaboratoryAngleEnergy
openmc.data.AngleDistribution
openmc.data.EnergyDistribution
openmc.data.ArbitraryTabulated
@ -392,6 +395,24 @@ Angle-Energy Distributions
openmc.data.LevelInelastic
openmc.data.ContinuousTabular
Resonance Data
--------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.data.Resonances
openmc.data.ResonanceRange
openmc.data.SingleLevelBreitWigner
openmc.data.MultiLevelBreitWigner
openmc.data.ReichMoore
openmc.data.RMatrixLimited
openmc.data.ParticlePair
openmc.data.SpinGroup
openmc.data.Unresolved
ACE Format
----------
@ -412,9 +433,37 @@ Functions
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.data.ace.ascii_to_binary
openmc.data.write_compact_458_library
ENDF Format
-----------
Classes
+++++++
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.data.endf.Evaluation
Functions
+++++++++
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.data.endf.float_endf
openmc.data.endf.get_cont_record
openmc.data.endf.get_head_record
openmc.data.endf.get_tab1_record
openmc.data.endf.get_tab2_record
openmc.data.endf.get_text_record
.. _Jupyter: https://jupyter.org/
.. _NumPy: http://www.numpy.org/

2
openmc/data/__init__.py
View file

@ -4,6 +4,7 @@ from .reaction import *
from .ace import *
from .angle_distribution import *
from .function import *
from .endf import *
from .energy_distribution import *
from .product import *
from .angle_energy import *
@ -15,3 +16,4 @@ from .thermal import *
from .urr import *
from .library import *
from .fission_energy import *
from .resonance import *

99
openmc/data/angle_distribution.py
View file

@ -1,12 +1,15 @@
from collections import Iterable
from io import StringIO
from numbers import Real
import numpy as np
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Univariate, Tabular, Uniform
from openmc.stats import Univariate, Tabular, Uniform, Legendre
from .function import INTERPOLATION_SCHEME
from .endf import get_head_record, get_cont_record, get_tab1_record, \
get_list_record, get_tab2_record
class AngleDistribution(EqualityMixin):
@ -199,3 +202,97 @@ class AngleDistribution(EqualityMixin):
mu.append(mu_i)
return cls(energy, mu)
@classmethod
def from_endf(cls, ev, mt):
"""Generate an angular distribution from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
mt : int
The MT value of the reaction to get angular distributions for
Returns
-------
openmc.data.AngleDistribution
Angular distribution
"""
file_obj = StringIO(ev.section[4, mt])
# Read HEAD record
items = get_head_record(file_obj)
ltt = items[3]
# Read CONT record
items = get_cont_record(file_obj)
li = items[2]
center_of_mass = (items[3] == 2)
if ltt == 0 and li == 1:
# Purely isotropic
energy = np.array([0., ev.info['energy_max']])
mu = [Uniform(-1., 1.), Uniform(-1., 1.)]
elif ltt == 1 and li == 0:
# Legendre polynomial coefficients
params, tab2 = get_tab2_record(file_obj)
n_energy = params[5]
energy = np.zeros(n_energy)
mu = []
for i in range(n_energy):
items, al = get_list_record(file_obj)
temperature = items[0]
energy[i] = items[1]
coefficients = np.asarray([1.0] + al)
mu.append(Legendre(coefficients))
elif ltt == 2 and li == 0:
# Tabulated probability distribution
params, tab2 = get_tab2_record(file_obj)
n_energy = params[5]
energy = np.zeros(n_energy)
mu = []
for i in range(n_energy):
params, f = get_tab1_record(file_obj)
temperature = params[0]
energy[i] = params[1]
if f.n_regions > 1:
raise NotImplementedError('Angular distribution with multiple '
'interpolation regions not supported.')
mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))
elif ltt == 3 and li == 0:
# Legendre for low energies / tabulated for high energies
params, tab2 = get_tab2_record(file_obj)
n_energy_legendre = params[5]
energy_legendre = np.zeros(n_energy_legendre)
mu = []
for i in range(n_energy_legendre):
items, al = get_list_record(file_obj)
temperature = items[0]
energy_legendre[i] = items[1]
coefficients = np.asarray([1.0] + al)
mu.append(Legendre(coefficients))
params, tab2 = get_tab2_record(file_obj)
n_energy_tabulated = params[5]
energy_tabulated = np.zeros(n_energy_tabulated)
for i in range(n_energy_tabulated):
params, f = get_tab1_record(file_obj)
temperature = params[0]
energy_tabulated[i] = params[1]
if f.n_regions > 1:
raise NotImplementedError('Angular distribution with multiple '
'interpolation regions not supported.')
mu.append(Tabular(f.x, f.y, INTERPOLATION_SCHEME[f.interpolation[0]]))
energy = np.concatenate((energy_legendre, energy_tabulated))
return AngleDistribution(energy, mu)

3
openmc/data/angle_energy.py
View file

@ -1,4 +1,5 @@
from abc import ABCMeta, abstractmethod
from io import StringIO
import openmc.data
from openmc.mixin import EqualityMixin
@ -40,7 +41,7 @@ class AngleEnergy(EqualityMixin):
@staticmethod
def from_ace(ace, location_dist, location_start, rx=None):
"""Generate an AngleEnergy object from ACE data
"""Generate an angle-energy distribution from ACE data
Parameters
----------

53
openmc/data/correlated.py
View file

@ -5,9 +5,11 @@ from warnings import warn
import numpy as np
import openmc.checkvalue as cv
from openmc.stats import Tabular, Univariate, Discrete, Mixture, Uniform
from openmc.stats import Tabular, Univariate, Discrete, Mixture, \
Uniform, Legendre
from .function import INTERPOLATION_SCHEME
from .angle_energy import AngleEnergy
from .endf import get_list_record, get_tab2_record
class CorrelatedAngleEnergy(AngleEnergy):
@ -405,3 +407,52 @@ class CorrelatedAngleEnergy(AngleEnergy):
mu.append(mu_i)
return cls(breakpoints, interpolation, energy, energy_out, mu)
@classmethod
def from_endf(cls, file_obj):
"""Generate correlated angle-energy distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for a correlated
angle-energy distribution
Returns
-------
openmc.data.CorrelatedAngleEnergy
Correlated angle-energy distribution
"""
params, tab2 = get_tab2_record(file_obj)
lep = params[3]
ne = params[5]
energy = np.zeros(ne)
n_discrete_energies = np.zeros(ne, dtype=int)
energy_out = []
mu = []
for i in range(ne):
items, values = get_list_record(file_obj)
energy[i] = items[1]
n_discrete_energies[i] = items[2]
# TODO: separate out discrete lines
n_angle = items[3]
n_energy_out = items[5]
values = np.asarray(values)
values.shape = (n_energy_out, n_angle + 2)
# Outgoing energy distribution at the i-th incoming energy
eout_i = values[:,0]
eout_p_i = values[:,1]
energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep],
ignore_negative=True)
energy_out.append(energy_out_i)
# Legendre coefficients used for angular distributions
mu_i = []
for j in range(n_energy_out):
mu_i.append(Legendre(values[j,1:]))
mu.append(mu_i)
return cls(tab2.breakpoints, tab2.interpolation, energy,
energy_out, mu)

60
openmc/data/data.py
View file

@ -104,8 +104,8 @@ NATURAL_ABUNDANCE = {
'U234': 5.4e-05, 'U235': 0.007204, 'U238': 0.992742
}
ATOMIC_SYMBOL = {1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N',
8: 'O', 9: 'F', 10: 'Ne', 11: 'Na', 12: 'Mg', 13: 'Al',
ATOMIC_SYMBOL = {0: 'n', 1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C',
7: 'N', 8: 'O', 9: 'F', 10: 'Ne', 11: 'Na', 12: 'Mg', 13: 'Al',
14: 'Si', 15: 'P', 16: 'S', 17: 'Cl', 18: 'Ar', 19: 'K',
20: 'Ca', 21: 'Sc', 22: 'Ti', 23: 'V', 24: 'Cr', 25: 'Mn',
26: 'Fe', 27: 'Co', 28: 'Ni', 29: 'Cu', 30: 'Zn', 31: 'Ga',
@ -129,60 +129,6 @@ ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
_ATOMIC_MASS = {}
REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
5: '(n,misc)', 11: '(n,2nd)', 16: '(n,2n)', 17: '(n,3n)',
18: '(n,fission)', 19: '(n,f)', 20: '(n,nf)', 21: '(n,2nf)',
22: '(n,na)', 23: '(n,n3a)', 24: '(n,2na)', 25: '(n,3na)',
27: '(n,absorption)', 28: '(n,np)', 29: '(n,n2a)',
30: '(n,2n2a)', 32: '(n,nd)', 33: '(n,nt)', 34: '(n,nHe-3)',
35: '(n,nd2a)', 36: '(n,nt2a)', 37: '(n,4n)', 38: '(n,3nf)',
41: '(n,2np)', 42: '(n,3np)', 44: '(n,n2p)', 45: '(n,npa)',
91: '(n,nc)', 101: '(n,disappear)', 102: '(n,gamma)',
103: '(n,p)', 104: '(n,d)', 105: '(n,t)', 106: '(n,3He)',
107: '(n,a)', 108: '(n,2a)', 109: '(n,3a)', 111: '(n,2p)',
112: '(n,pa)', 113: '(n,t2a)', 114: '(n,d2a)', 115: '(n,pd)',
116: '(n,pt)', 117: '(n,da)', 152: '(n,5n)', 153: '(n,6n)',
154: '(n,2nt)', 155: '(n,ta)', 156: '(n,4np)', 157: '(n,3nd)',
158: '(n,nda)', 159: '(n,2npa)', 160: '(n,7n)', 161: '(n,8n)',
162: '(n,5np)', 163: '(n,6np)', 164: '(n,7np)', 165: '(n,4na)',
166: '(n,5na)', 167: '(n,6na)', 168: '(n,7na)', 169: '(n,4nd)',
170: '(n,5nd)', 171: '(n,6nd)', 172: '(n,3nt)', 173: '(n,4nt)',
174: '(n,5nt)', 175: '(n,6nt)', 176: '(n,2n3He)',
177: '(n,3n3He)', 178: '(n,4n3He)', 179: '(n,3n2p)',
180: '(n,3n3a)', 181: '(n,3npa)', 182: '(n,dt)',
183: '(n,npd)', 184: '(n,npt)', 185: '(n,ndt)',
186: '(n,np3He)', 187: '(n,nd3He)', 188: '(n,nt3He)',
189: '(n,nta)', 190: '(n,2n2p)', 191: '(n,p3He)',
192: '(n,d3He)', 193: '(n,3Hea)', 194: '(n,4n2p)',
195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)'}
REACTION_NAME.update({i: '(n,n{})'.format(i-50) for i in range(50, 91)})
REACTION_NAME.update({i: '(n,p{})'.format(i-600) for i in range(600, 649)})
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))}
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
@ -207,7 +153,7 @@ def atomic_mass(isotope):
mass_file = os.path.join(os.path.dirname(__file__), 'mass.mas12')
with open(mass_file, 'r') as ame:
# Read lines in file starting at line 40
for line in itertools.islice(ame, 40, None):
for line in itertools.islice(ame, 39, None):
name = '{}{}'.format(line[20:22].strip(), int(line[16:19]))
mass = float(line[96:99]) + 1e-6*float(
line[100:106] + '.' + line[107:112])

409
openmc/data/endf.py Normal file
View file

@ -0,0 +1,409 @@
"""Module for parsing and manipulating data from ENDF evaluations.
All the classes and functions in this module are based on document
ENDF-102 titled "Data Formats and Procedures for the Evaluated Nuclear
Data File ENDF-6". The latest version from June 2009 can be found at
http://www-nds.iaea.org/ndspub/documents/endf/endf102/endf102.pdf
"""
from __future__ import print_function, division, unicode_literals
import io
import re
import os
from math import pi
from collections import OrderedDict, Iterable
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from .function import Tabulated1D, INTERPOLATION_SCHEME
from openmc.stats.univariate import Uniform, Tabular, Legendre
LIBRARIES = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF',
4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL',
31: 'INDL/V', 32: 'INDL/A', 33: 'FENDL', 34: 'IRDF',
35: 'BROND', 36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'}
SUM_RULES = {1: [2, 3],
3: [4, 5, 11, 16, 17, 22, 23, 24, 25, 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: [8, 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))}
_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
def radiation_type(value):
p = {0: 'gamma', 1: 'beta-', 2: 'ec/beta+', 3: 'IT',
4: 'alpha', 5: 'neutron', 6: 'sf', 7: 'proton',
8: 'e-', 9: 'xray', 10: 'unknown'}
if value % 1.0 == 0:
return p[int(value)]
else:
return (p[int(value)], p[int(10*value % 10)])
def float_endf(s):
"""Convert string of floating point number in ENDF to float.
The ENDF-6 format uses an 'e-less' floating point number format,
e.g. -1.23481+10. Trying to convert using the float built-in won't work
because of the lack of an 'e'. This function allows such strings to be
converted while still allowing numbers that are not in exponential notation
to be converted as well.
Parameters
----------
s : str
Floating-point number from an ENDF file
Returns
-------
float
The number
"""
return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s))
def get_text_record(file_obj):
"""Return data from a TEXT record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
str
Text within the TEXT record
"""
return file_obj.readline()[:66]
def get_cont_record(file_obj, skipC=False):
"""Return data from a CONT record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
list
The six items within the CONT record
"""
line = file_obj.readline()
if skipC:
C1 = None
C2 = None
else:
C1 = float_endf(line[:11])
C2 = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
N1 = int(line[44:55])
N2 = int(line[55:66])
return [C1, C2, L1, L2, N1, N2]
def get_head_record(file_obj):
"""Return data from a HEAD record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
list
The six items within the HEAD record
"""
line = file_obj.readline()
ZA = int(float_endf(line[:11]))
AWR = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
N1 = int(line[44:55])
N2 = int(line[55:66])
return [ZA, AWR, L1, L2, N1, N2]
def get_list_record(file_obj):
"""Return data from a LIST record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
list
The six items within the header
list
The values within the list
"""
# determine how many items are in list
items = get_cont_record(file_obj)
NPL = items[4]
# read items
b = []
for i in range((NPL - 1)//6 + 1):
line = file_obj.readline()
n = min(6, NPL - 6*i)
for j in range(n):
b.append(float_endf(line[11*j:11*(j + 1)]))
return (items, b)
def get_tab1_record(file_obj):
"""Return data from a TAB1 record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
Returns
-------
list
The six items within the header
openmc.data.Tabulated1D
The tabulated function
"""
# Determine how many interpolation regions and total points there are
line = file_obj.readline()
C1 = float_endf(line[:11])
C2 = float_endf(line[11:22])
L1 = int(line[22:33])
L2 = int(line[33:44])
n_regions = int(line[44:55])
n_pairs = int(line[55:66])
params = [C1, C2, L1, L2]
# Read the interpolation region data, namely NBT and INT
breakpoints = np.zeros(n_regions, dtype=int)
interpolation = np.zeros(n_regions, dtype=int)
m = 0
for i in range((n_regions - 1)//3 + 1):
line = file_obj.readline()
to_read = min(3, n_regions - m)
for j in range(to_read):
breakpoints[m] = int(line[0:11])
interpolation[m] = int(line[11:22])
line = line[22:]
m += 1
# Read tabulated pairs x(n) and y(n)
x = np.zeros(n_pairs)
y = np.zeros(n_pairs)
m = 0
for i in range((n_pairs - 1)//3 + 1):
line = file_obj.readline()
to_read = min(3, n_pairs - m)
for j in range(to_read):
x[m] = float_endf(line[:11])
y[m] = float_endf(line[11:22])
line = line[22:]
m += 1
return params, Tabulated1D(x, y, breakpoints, interpolation)
def get_tab2_record(file_obj):
# Determine how many interpolation regions and total points there are
params = get_cont_record(file_obj)
n_regions = params[4]
# Read the interpolation region data, namely NBT and INT
breakpoints = np.zeros(n_regions, dtype=int)
interpolation = np.zeros(n_regions, dtype=int)
m = 0
for i in range((n_regions - 1)//3 + 1):
line = file_obj.readline()
to_read = min(3, n_regions - m)
for j in range(to_read):
breakpoints[m] = int(line[0:11])
interpolation[m] = int(line[11:22])
line = line[22:]
m += 1
return params, Tabulated2D(breakpoints, interpolation)
class Evaluation(object):
"""ENDF material evaluation with multiple files/sections
Parameters
----------
filename : str
Path to ENDF file to read
Attributes
----------
info : dict
Miscallaneous information about the evaluation.
target : dict
Information about the target material, such as its mass, isomeric state,
whether it's stable, and whether it's fissionable.
projectile : dict
Information about the projectile such as its mass.
reaction_list : list of 4-tuples
List of sections in the evaluation. The entries of the tuples are the
file (MF), section (MT), number of records (NC), and modification
indicator (MOD).
"""
def __init__(self, filename):
fh = open(filename, 'r')
self.section = {}
self.info = {}
self.target = {}
self.projectile = {}
self.reaction_list = []
# Determine MAT number for this evaluation
MF = 0
while MF == 0:
position = fh.tell()
line = fh.readline()
MF = int(line[70:72])
self.material = int(line[66:70])
fh.seek(position)
while True:
# Find next section
while True:
position = fh.tell()
line = fh.readline()
MAT = int(line[66:70])
MF = int(line[70:72])
MT = int(line[72:75])
if MT > 0 or MAT == 0:
fh.seek(position)
break
# If end of material reached, exit loop
if MAT == 0:
break
section_data = ''
while True:
line = fh.readline()
if line[72:75] == ' 0':
break
else:
section_data += line
self.section[MF, MT] = section_data
self._read_header()
def _read_header(self):
file_obj = io.StringIO(self.section[1, 451])
# Information about target/projectile
items = get_head_record(file_obj)
self.target['atomic_number'] = items[0] // 1000
self.target['mass_number'] = items[0] % 1000
self.target['mass'] = items[1]
self._LRP = items[2]
self.target['fissionable'] = (items[3] == 1)
try:
global LIBRARIES
library = LIBRARIES[items[4]]
except KeyError:
library = 'Unknown'
self.info['modification'] = items[5]
# Control record 1
items = get_cont_record(file_obj)
self.target['excitation_energy'] = items[0]
self.target['stable'] = (int(items[1]) == 0)
self.target['state'] = items[2]
self.target['isomeric_state'] = items[3]
self.info['format'] = items[5]
assert self.info['format'] == 6
# Control record 2
items = get_cont_record(file_obj)
self.projectile['mass'] = items[0]
self.info['energy_max'] = items[1]
library_release = items[2]
self.info['sublibrary'] = items[4]
library_version = items[5]
self.info['library'] = (library, library_version, library_release)
# Control record 3
items = get_cont_record(file_obj)
self.target['temperature'] = items[0]
self.info['derived'] = (items[2] > 0)
NWD = items[4]
NXC = items[5]
# Text records
text = [get_text_record(file_obj) for i in range(NWD)]
if len(text) >= 5:
self.target['zsymam'] = text[0][0:11]
self.info['laboratory'] = text[0][11:22]
self.info['date'] = text[0][22:32]
self.info['author'] = text[0][32:66]
self.info['reference'] = text[1][1:22]
self.info['date_distribution'] = text[1][22:32]
self.info['date_release'] = text[1][33:43]
self.info['date_entry'] = text[1][55:63]
self.info['identifier'] = text[2:5]
self.info['description'] = text[5:]
# File numbers, reaction designations, and number of records
for i in range(NXC):
items = get_cont_record(file_obj, skipC=True)
MF, MT, NC, MOD = items[2:6]
self.reaction_list.append((MF, MT, NC, MOD))
class Tabulated2D(object):
"""Metadata for a two-dimensional function.
This is a dummy class that is not really used other than to store the
interpolation information for a two-dimensional function. Once we refactor
to adopt GND-like data containers, this will probably be removed or
extended.
Parameters
----------
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).
"""
def __init__(self, breakpoints, interpolation):
self.breakpoints = breakpoints
self.interpolation = interpolation

44
openmc/data/endf_utils.py
View file

@ -1,44 +0,0 @@
"""This module contains a few utility functions for reading ENDF_ data. It is by
no means enough to read an entire ENDF file. For a more complete ENDF reader,
see Pyne_.
.. _ENDF: http://www.nndc.bnl.gov/endf
.. _Pyne: http://www.pyne.io
"""
import re
def read_float(float_string):
"""Parse ENDF 6E11.0 formatted string into a float."""
assert len(float_string) == 11
pattern = r'([\s\-]\d+\.\d+)([\+\-]\d+)'
return float(re.sub(pattern, r'\1e\2', float_string))
def read_CONT_line(line):
"""Parse 80-column line from ENDF CONT record into floats and ints."""
return (read_float(line[0:11]), read_float(line[11:22]), int(line[22:33]),
int(line[33:44]), int(line[44:55]), int(line[55:66]),
int(line[66:70]), int(line[70:72]), int(line[72:75]),
int(line[75:80]))
def identify_nuclide(fname):
"""Read the header of an ENDF file and extract identifying information."""
with open(fname, 'r') as fh:
# Skip the tape id (TPID).
line = fh.readline()
# Read the first HEAD and CONT info.
line = fh.readline()
ZA, AW, LRP, LFI, NLIB, NMOD, MAT, MF, MT, NS = read_CONT_line(line)
line = fh.readline()
ELIS, STA, LIS, LISO, junk, NFOR, MAT, MF, MT, NS = read_CONT_line(line)
# Return dictionary of the most important identifying information.
return {'Z': int(ZA) // 1000,
'A': int(ZA) % 1000,
'LFI': bool(LFI),
'LIS': LIS,
'LISO': LISO}

189
openmc/data/energy_distribution.py
View file

@ -9,6 +9,7 @@ from .function import Tabulated1D, INTERPOLATION_SCHEME
from openmc.stats.univariate import Univariate, Tabular, Discrete, Mixture
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .endf import get_tab1_record, get_tab2_record
class EnergyDistribution(EqualityMixin):
@ -57,6 +58,40 @@ class EnergyDistribution(EqualityMixin):
raise ValueError("Unknown energy distribution type: {}"
.format(energy_type))
@staticmethod
def from_endf(file_obj, params):
"""Generate energy distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.EnergyDistribution
A sub-class of :class:`openmc.data.EnergyDistribution`
"""
lf = params[3]
if lf == 1:
return ArbitraryTabulated.from_endf(file_obj, params)
elif lf == 5:
return GeneralEvaporation.from_endf(file_obj, params)
elif lf == 7:
return MaxwellEnergy.from_endf(file_obj, params)
elif lf == 9:
return Evaporation.from_endf(file_obj, params)
elif lf == 11:
return WattEnergy.from_endf(file_obj, params)
elif lf == 12:
return MadlandNix.from_endf(file_obj, params)
class ArbitraryTabulated(EnergyDistribution):
r"""Arbitrary tabulated function given in ENDF MF=5, LF=1 represented as
@ -88,6 +123,37 @@ class ArbitraryTabulated(EnergyDistribution):
def to_hdf5(self, group):
raise NotImplementedError
@classmethod
def from_endf(cls, file_obj, params):
"""Generate arbitrary tabulated distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.ArbitraryTabulated
Arbitrary tabulated distribution
"""
params, tab2 = get_tab2_record(file_obj)
n_energies = params[5]
energy = np.zeros(n_energies)
pdf = []
for j in range(n_energies):
params, func = get_tab1_record(file_obj)
energy[j] = params[1]
pdf.append(func)
return cls(energy, pdf)
class GeneralEvaporation(EnergyDistribution):
r"""General evaporation spectrum given in ENDF MF=5, LF=5 represented as
@ -130,6 +196,31 @@ class GeneralEvaporation(EnergyDistribution):
def from_ace(cls, ace, idx=0):
raise NotImplementedError
@classmethod
def from_endf(cls, file_obj, params):
"""Generate general evaporation spectrum from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.GeneralEvaporation
General evaporation spectrum
"""
u = params[0]
params, theta = get_tab1_record(file_obj)
params, g = get_tab1_record(file_obj)
return cls(theta, g, u)
class MaxwellEnergy(EnergyDistribution):
r"""Simple Maxwellian fission spectrum represented as
@ -238,6 +329,30 @@ class MaxwellEnergy(EnergyDistribution):
return cls(theta, u)
@classmethod
def from_endf(cls, file_obj, params):
"""Generate Maxwell distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.MaxwellEnergy
Maxwell distribution
"""
u = params[0]
params, theta = get_tab1_record(file_obj)
return cls(theta, u)
class Evaporation(EnergyDistribution):
r"""Evaporation spectrum represented as
@ -346,6 +461,30 @@ class Evaporation(EnergyDistribution):
return cls(theta, u)
@classmethod
def from_endf(cls, file_obj, params):
"""Generate evaporation spectrum from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.Evaporation
Evaporation spectrum
"""
u = params[0]
params, theta = get_tab1_record(file_obj)
return cls(theta, u)
class WattEnergy(EnergyDistribution):
r"""Energy-dependent Watt spectrum represented as
@ -480,6 +619,31 @@ class WattEnergy(EnergyDistribution):
return cls(a, b, u)
@classmethod
def from_endf(cls, file_obj, params):
"""Generate Watt fission spectrum from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.WattEnergy
Watt fission spectrum
"""
u = params[0]
params, a = get_tab1_record(file_obj)
params, b = get_tab1_record(file_obj)
return cls(a, b, u)
class MadlandNix(EnergyDistribution):
r"""Energy-dependent fission neutron spectrum (Madland and Nix) given in
@ -587,6 +751,31 @@ class MadlandNix(EnergyDistribution):
tm = Tabulated1D.from_hdf5(group['tm'])
return cls(efl, efh, tm)
@classmethod
def from_endf(cls, file_obj, params):
"""Generate Madland-Nix fission spectrum from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for an energy
distribution.
params : list
List of parameters at the start of the energy distribution that
includes the LF value indicating what type of energy distribution is
present.
Returns
-------
openmc.data.MadlandNix
Madland-Nix fission spectrum
"""
params, tm = get_tab1_record(file_obj)
efl, efh = params[0:2]
return cls(efl, efh, tm)
class DiscretePhoton(EnergyDistribution):
"""Discrete photon energy distribution

68
openmc/data/fission_energy.py
View file

@ -1,12 +1,13 @@
from collections import Callable
from copy import deepcopy
from io import StringIO
import sys
import h5py
import numpy as np
from .data import ATOMIC_SYMBOL
from .endf_utils import read_float, read_CONT_line, identify_nuclide
from .endf import get_cont_record, get_list_record, Evaluation
from .function import Function1D, Tabulated1D, Polynomial, Sum
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -15,13 +16,13 @@ if sys.version_info[0] >= 3:
basestring = str
def _extract_458_data(filename):
def _extract_458_data(ev):
"""Read an ENDF file and extract the MF=1, MT=458 values.
Parameters
----------
filename : str
Path to and ENDF file
ev : openmc.data.Evaluation
ENDF evaluation
Returns
-------
@ -37,33 +38,22 @@ def _extract_458_data(filename):
caution.
"""
ident = identify_nuclide(filename)
if not ident['LFI']:
if not ev.target['fissionable']:
# This nuclide isn't fissionable.
return None
# Extract the MF=1, MT=458 section.
lines = []
with open(filename, 'r') as fh:
line = fh.readline()
while line != '':
if line[70:75] == ' 1458':
lines.append(line)
line = fh.readline()
if len(lines) == 0:
if (1, 458) not in ev.section:
# No 458 data here.
return None
file_obj = StringIO(ev.section[1, 458])
# Read the number of coefficients in this LIST record.
NPL = read_CONT_line(lines[1])[4]
items = get_cont_record(file_obj)
NPL = items[3]
# Parse the ENDF LIST into an array.
data = []
for i in range(NPL):
row, column = divmod(i, 6)
data.append(read_float(lines[2 + row][11*column:11*(column+1)]))
items, data = get_list_record(file_obj)
# Declare the coefficient names and the order they are given in. The LIST
# contains a value followed immediately by an uncertainty for each of these
@ -161,20 +151,22 @@ def write_compact_458_library(endf_files, output_name='fission_Q_data.h5',
for fname in endf_files:
if verbose: print(fname)
ident = identify_nuclide(fname)
ev = Evaluation(fname)
# Skip non-fissionable nuclides.
if not ident['LFI']: continue
if not ev.target['fissionable']:
continue
# Get the important bits.
data = _extract_458_data(fname)
data = _extract_458_data(ev)
if data is None: continue
value, uncertainty = data
# Make a group for this isomer.
name = ATOMIC_SYMBOL[ident['Z']] + str(ident['A'])
if ident['LISO'] != 0:
name += '_m' + str(ident['LISO'])
name = ATOMIC_SYMBOL[ev.target['atomic_number']] + \
str(ev.target['mass_number'])
if ev.target['isomeric_state'] != 0:
name += '_m' + str(ev.target['isomeric_state'])
nuclide_group = out.create_group(name)
# Write all the coefficients into one array. The first dimension gives
@ -371,7 +363,6 @@ class FissionEnergyRelease(EqualityMixin):
component. The keys are the 2-3 letter strings used in ENDF-102,
e.g. 'EFR' and 'ET'. The list will have a length of 1 for Sher-Beck
data, more for polynomial data.
incident_neutron : openmc.data.IncidentNeutron
Corresponding incident neutron dataset
@ -440,14 +431,13 @@ class FissionEnergyRelease(EqualityMixin):
return out
@classmethod
def from_endf(cls, filename, incident_neutron):
def from_endf(cls, ev, incident_neutron):
"""Generate fission energy release data from an ENDF file.
Parameters
----------
filename : str
Name of the ENDF file containing fission energy release data
ev : openmc.data.endf.Evaluation
ENDF evaluation
incident_neutron : openmc.data.IncidentNeutron
Corresponding incident neutron dataset
@ -459,21 +449,20 @@ class FissionEnergyRelease(EqualityMixin):
"""
# Check to make sure this ENDF file matches the expected isomer.
ident = identify_nuclide(filename)
if ident['Z'] != incident_neutron.atomic_number:
if ev.target['atomic_number'] != incident_neutron.atomic_number:
raise ValueError('The atomic number of the ENDF evaluation does '
'not match the given IncidentNeutron.')
if ident['A'] != incident_neutron.mass_number:
if ev.target['mass_number'] != incident_neutron.mass_number:
raise ValueError('The atomic mass of the ENDF evaluation does '
'not match the given IncidentNeutron.')
if ident['LISO'] != incident_neutron.metastable:
if ev.target['isomeric_state'] != incident_neutron.metastable:
raise ValueError('The metastable state of the ENDF evaluation does '
'not match the given IncidentNeutron.')
if not ident['LFI']:
if not ev.target['fissionable']:
raise ValueError('The ENDF evaluation is not fissionable.')
# Read the 458 data from the ENDF file.
value, uncertainty = _extract_458_data(filename)
value, uncertainty = _extract_458_data(ev)
# Build the object.
return cls._from_dictionary(value, incident_neutron)
@ -516,7 +505,6 @@ class FissionEnergyRelease(EqualityMixin):
Path to an HDF5 file containing fission energy release data. This
file should have been generated form the
:func:`openmc.data.write_compact_458_library` function.
incident_neutron : openmc.data.IncidentNeutron
Corresponding incident neutron dataset

78
openmc/data/function.py
View file

@ -4,6 +4,7 @@ from numbers import Real, Integral
import numpy as np
import openmc.data
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -423,3 +424,80 @@ class Sum(EqualityMixin):
def functions(self, functions):
cv.check_type('functions', functions, Iterable, Callable)
self._functions = functions
class ResonancesWithBackground(EqualityMixin):
"""Cross section in resolved resonance region.
Parameters
----------
resonances : openmc.data.Resonances
Resolved resonance parameter data
background : Callable
Background cross section as a function of energy
mt : int
MT value of the reaction
Attributes
----------
resonances : openmc.data.Resonances
Resolved resonance parameter data
background : Callable
Background cross section as a function of energy
mt : int
MT value of the reaction
"""
def __init__(self, resonances, background, mt):
self.resonances = resonances
self.background = background
self.mt = mt
def __call__(self, x):
# Get background cross section
xs = self.background(x)
rrr = self.resonances.resolved
if isinstance(x, Iterable):
# Determine which energies are within resolved resonance range
within_rrr = (rrr.energy_min <= x) & (x <= rrr.energy_max)
# Get resonance cross sections and add to background
resonant_xs = rrr.reconstruct(x[within_rrr])
xs[within_rrr] += resonant_xs[self.mt]
else:
if rrr.energy_min <= x <= rrr.energy_max:
resonant_xs = rrr.reconstruct(x)
xs += resonant_xs[self.mt]
return xs
@property
def background(self):
return self._background
@property
def mt(self):
return self._mt
@property
def resonances(self):
return self._resonances
@background.setter
def background(self, background):
cv.check_type('background cross section', background, Callable)
self._background = background
@mt.setter
def mt(self, mt):
cv.check_type('MT value', mt, Integral)
self._mt = mt
@resonances.setter
def resonances(self, resonances):
cv.check_type('resolved resonance parameters', resonances,
openmc.data.Resonances)
self._resonances = resonances

52
openmc/data/kalbach_mann.py
View file

@ -8,6 +8,7 @@ import openmc.checkvalue as cv
from openmc.stats import Tabular, Univariate, Discrete, Mixture
from .function import Tabulated1D, INTERPOLATION_SCHEME
from .angle_energy import AngleEnergy
from .endf import get_list_record, get_tab2_record
class KalbachMann(AngleEnergy):
@ -346,3 +347,54 @@ class KalbachMann(AngleEnergy):
km_a.append(Tabulated1D(data[0], data[4]))
return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
@classmethod
def from_endf(cls, file_obj):
"""Generate Kalbach-Mann distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of the Kalbach-Mann distribution
Returns
-------
openmc.data.KalbachMann
Kalbach-Mann energy-angle distribution
"""
params, tab2 = get_tab2_record(file_obj)
lep = params[3]
ne = params[5]
energy = np.zeros(ne)
n_discrete_energies = np.zeros(ne, dtype=int)
energy_out = []
precompound = []
slope = []
for i in range(ne):
items, values = get_list_record(file_obj)
energy[i] = items[1]
n_discrete_energies[i] = items[2]
# TODO: split out discrete energies
n_angle = items[3]
n_energy_out = items[5]
values = np.asarray(values)
values.shape = (n_energy_out, n_angle + 2)
# Outgoing energy distribution at the i-th incoming energy
eout_i = values[:,0]
eout_p_i = values[:,1]
energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep])
energy_out.append(energy_out_i)
# Precompound and slope factors for Kalbach-Mann
r_i = values[:,2]
if n_angle == 2:
a_i = values[:,3]
else:
a_i = np.zeros_like(r_i)
precompound.append(Tabulated1D(eout_i, r_i))
slope.append(Tabulated1D(eout_i, a_i))
return cls(tab2.breakpoints, tab2.interpolation, energy,
energy_out, precompound, slope)

139
openmc/data/laboratory.py Normal file
View file

@ -0,0 +1,139 @@
from collections import Iterable
from numbers import Real, Integral
import numpy as np
import openmc.checkvalue as cv
from openmc.stats import Tabular, Univariate, Discrete, Mixture
from .angle_energy import AngleEnergy
from .function import INTERPOLATION_SCHEME
from .endf import get_tab2_record, get_tab1_record
class LaboratoryAngleEnergy(AngleEnergy):
"""Laboratory angle-energy distribution
Parameters
----------
breakpoints : Iterable of int
Breakpoints defining interpolation regions
interpolation : Iterable of int
Interpolation codes
energy : Iterable of float
Incoming energies at which distributions exist
mu : Iterable of openmc.stats.Univariate
Distribution of scattering cosines for each incoming energy
energy_out : Iterable of Iterable of openmc.stats.Univariate
Distribution of outgoing energies for each incoming energy/scattering
cosine
Attributes
----------
breakpoints : Iterable of int
Breakpoints defining interpolation regions
interpolation : Iterable of int
Interpolation codes
energy : Iterable of float
Incoming energies at which distributions exist
mu : Iterable of openmc.stats.Univariate
Distribution of scattering cosines for each incoming energy
energy_out : Iterable of Iterable of openmc.stats.Univariate
Distribution of outgoing energies for each incoming energy/scattering
cosine
"""
def __init__(self, breakpoints, interpolation, energy, mu, energy_out):
super(LaboratoryAngleEnergy).__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.energy = energy
self.mu = mu
self.energy_out = energy_out
@property
def breakpoints(self):
return self._breakpoints
@property
def interpolation(self):
return self._interpolation
@property
def energy(self):
return self._energy
@property
def mu(self):
return self._mu
@property
def energy_out(self):
return self._energy_out
@breakpoints.setter
def breakpoints(self, breakpoints):
cv.check_type('laboratory angle-energy breakpoints', breakpoints,
Iterable, Integral)
self._breakpoints = breakpoints
@interpolation.setter
def interpolation(self, interpolation):
cv.check_type('laboratory angle-energy interpolation', interpolation,
Iterable, Integral)
self._interpolation = interpolation
@energy.setter
def energy(self, energy):
cv.check_type('laboratory angle-energy incoming energy', energy,
Iterable, Real)
self._energy = energy
@mu.setter
def mu(self, mu):
cv.check_type('laboratory angle-energy outgoing cosine', mu,
Iterable, Univariate)
self._mu = mu
@energy_out.setter
def energy_out(self, energy_out):
cv.check_iterable_type('laboratory angle-energy outgoing energy',
energy_out, Univariate, 2, 2)
self._energy_out = energy_out
@classmethod
def from_endf(cls, file_obj):
"""Generate laboratory angle-energy distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positioned at the start of a section for a correlated
angle-energy distribution
Returns
-------
openmc.data.LaboratoryAngleEnergy
Laboratory angle-energy distribution
"""
params, tab2 = get_tab2_record(file_obj)
ne = params[5]
energy = np.zeros(ne)
mu = []
energy_out = []
for i in range(ne):
params, tab2mu = get_tab2_record(file_obj)
energy[i] = params[1]
n_mu = params[5]
mu_i = np.zeros(n_mu)
p_mu_i = np.zeros(n_mu)
energy_out_i = []
for j in range(n_mu):
params, f = get_tab1_record(file_obj)
mu_i[j] = params[1]
p_mu_i[j] = sum(f.y)
energy_out_i.append(Tabular(f.x, f.y))
mu.append(Tabular(mu_i, p_mu_i))
energy_out.append(energy_out_i)
return cls(tab2.breakpoints, tab2.interpolation, energy, mu, energy_out)

23
openmc/data/nbody.py
View file

@ -4,6 +4,7 @@ import numpy as np
import openmc.checkvalue as cv
from .angle_energy import AngleEnergy
from .endf import get_cont_record
class NBodyPhaseSpace(AngleEnergy):
"""N-body phase space distribution
@ -140,3 +141,25 @@ class NBodyPhaseSpace(AngleEnergy):
n_particles = int(ace.xss[idx])
total_mass = ace.xss[idx + 1]
return cls(total_mass, n_particles, ace.atomic_weight_ratio, q_value)
@classmethod
def from_endf(cls, file_obj):
"""Generate N-body phase space distribution from an ENDF evaluation
Parameters
----------
file_obj : file-like object
ENDF file positions at the start of the N-body phase space
distribution
Returns
-------
openmc.data.NBodyPhaseSpace
N-body phase space distribution
"""
items = get_cont_record(file_obj)
total_mass = items[0]
n_particles = items[5]
# TODO: get awr and Q value
return cls(total_mass, n_particles, 1.0, 0.0)

100
openmc/data/neutron.py
View file

@ -8,12 +8,14 @@ from warnings import warn
import numpy as np
import h5py
from .data import ATOMIC_SYMBOL, SUM_RULES, K_BOLTZMANN
from .ace import Table, get_table
from .data import ATOMIC_SYMBOL, K_BOLTZMANN
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum
from .function import Tabulated1D, Sum, ResonancesWithBackground
from .endf import Evaluation, SUM_RULES
from .product import Product
from .reaction import Reaction, _get_photon_products
from .reaction import Reaction, _get_photon_products_ace
from .resonance import Resonances, _RESOLVED
from .urr import ProbabilityTables
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
@ -137,6 +139,8 @@ class IncidentNeutron(EqualityMixin):
Contains the cross sections, secondary angle and energy distributions,
and other associated data for each reaction. The keys are the MT values
and the values are Reaction objects.
resonances : openmc.data.Resonances or None
Resonance parameters
summed_reactions : collections.OrderedDict
Contains summed cross sections, e.g., the total cross section. The keys
are the MT values and the values are Reaction objects.
@ -166,6 +170,7 @@ class IncidentNeutron(EqualityMixin):
self.reactions = OrderedDict()
self.summed_reactions = OrderedDict()
self._urr = {}
self._resonances = None
def __contains__(self, mt):
return mt in self.reactions or mt in self.summed_reactions
@ -212,6 +217,10 @@ class IncidentNeutron(EqualityMixin):
def reactions(self):
return self._reactions
@property
def resonances(self):
return self._resonances
@property
def summed_reactions(self):
return self._summed_reactions
@ -236,7 +245,7 @@ class IncidentNeutron(EqualityMixin):
@atomic_number.setter
def atomic_number(self, atomic_number):
cv.check_type('atomic number', atomic_number, Integral)
cv.check_greater_than('atomic number', atomic_number, 0)
cv.check_greater_than('atomic number', atomic_number, 0, True)
self._atomic_number = atomic_number
@mass_number.setter
@ -268,6 +277,11 @@ class IncidentNeutron(EqualityMixin):
cv.check_type('reactions', reactions, Mapping)
self._reactions = reactions
@resonances.setter
def resonances(self, resonances):
cv.check_type('resonances', resonances, Resonances)
self._resonances = resonances
@summed_reactions.setter
def summed_reactions(self, summed_reactions):
cv.check_type('summed reactions', summed_reactions, Mapping)
@ -587,7 +601,7 @@ class IncidentNeutron(EqualityMixin):
for mt_i in mts])
# Determine summed cross section
rx.products += _get_photon_products(ace, rx)
rx.products += _get_photon_products_ace(ace, rx)
data.summed_reactions[mt] = rx
# Read unresolved resonance probability tables
@ -596,3 +610,79 @@ class IncidentNeutron(EqualityMixin):
data.urr[strT] = urr
return data
@classmethod
def from_endf(cls, ev_or_filename):
"""Generate incident neutron continuous-energy data from an ENDF evaluation
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
Returns
-------
openmc.data.IncidentNeutron
Incident neutron continuous-energy data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
atomic_number = ev.target['atomic_number']
mass_number = ev.target['mass_number']
metastable = ev.target['isomeric_state']
atomic_weight_ratio = ev.target['mass']
temperature = ev.target['temperature']
# Determine name
element = ATOMIC_SYMBOL[atomic_number]
if metastable > 0:
name = '{}{}_m{}'.format(element, mass_number, metastable)
else:
name = '{}{}'.format(element, mass_number)
# Instantiate incident neutron data
data = cls(name, atomic_number, mass_number, metastable,
atomic_weight_ratio, temperature)
if (2, 151) in ev.section:
data.resonances = Resonances.from_endf(ev)
# Read each reaction
for mf, mt, nc, mod in ev.reaction_list:
if mf == 3:
data.reactions[mt] = Reaction.from_endf(ev, mt)
# Replace cross sections for elastic, capture, fission
try:
if isinstance(data.resonances.resolved, _RESOLVED):
for mt in (2, 102, 18):
if mt in data.reactions:
rx = data.reactions[mt]
rx.xs['0K'] = ResonancesWithBackground(
data.resonances, rx.xs['0K'], mt)
except ValueError:
# Thrown if multiple resolved ranges (e.g. Pu239 in ENDF/B-VII.1)
pass
# If first-chance, second-chance, etc. fission are present, check
# whether energy distributions were specified in MF=5. If not, copy the
# energy distribution from MT=18.
for mt, rx in data.reactions.items():
if mt in (19, 20, 21, 38):
if (5, mt) not in ev.section:
neutron = data.reactions[18].products[0]
rx.products[0].applicability = neutron.applicability
rx.products[0].distribution = neutron.distribution
# Read fission energy release (requires that we already know nu for
# fission)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)
data._evaluation = ev
return data

3
openmc/data/product.py
View file

@ -1,4 +1,5 @@
from collections import Iterable
from io import StringIO
from numbers import Real
import sys
@ -6,8 +7,8 @@ import numpy as np
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .function import Tabulated1D, Polynomial, Function1D
from .angle_energy import AngleEnergy
from .function import Tabulated1D, Polynomial, Function1D
if sys.version_info[0] >= 3:
basestring = str

493
openmc/data/reaction.py
View file

@ -3,21 +3,190 @@ from collections import Iterable, Callable, MutableMapping
from copy import deepcopy
from numbers import Real, Integral
from warnings import warn
from io import StringIO
import numpy as np
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Uniform
from openmc.stats import Uniform, Tabular, Legendre
from .angle_distribution import AngleDistribution
from .angle_energy import AngleEnergy
from .function import Tabulated1D, Polynomial, Function1D
from .data import REACTION_NAME, K_BOLTZMANN
from .correlated import CorrelatedAngleEnergy
from .data import ATOMIC_SYMBOL, K_BOLTZMANN
from .endf import get_head_record, get_tab1_record, get_list_record, \
get_tab2_record, get_cont_record
from .energy_distribution import EnergyDistribution, LevelInelastic, \
DiscretePhoton
from .function import Tabulated1D, Polynomial
from .kalbach_mann import KalbachMann
from .laboratory import LaboratoryAngleEnergy
from .nbody import NBodyPhaseSpace
from .product import Product
from .uncorrelated import UncorrelatedAngleEnergy
def _get_fission_products(ace):
REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
5: '(n,misc)', 11: '(n,2nd)', 16: '(n,2n)', 17: '(n,3n)',
18: '(n,fission)', 19: '(n,f)', 20: '(n,nf)', 21: '(n,2nf)',
22: '(n,na)', 23: '(n,n3a)', 24: '(n,2na)', 25: '(n,3na)',
27: '(n,absorption)', 28: '(n,np)', 29: '(n,n2a)',
30: '(n,2n2a)', 32: '(n,nd)', 33: '(n,nt)', 34: '(n,nHe-3)',
35: '(n,nd2a)', 36: '(n,nt2a)', 37: '(n,4n)', 38: '(n,3nf)',
41: '(n,2np)', 42: '(n,3np)', 44: '(n,n2p)', 45: '(n,npa)',
91: '(n,nc)', 101: '(n,disappear)', 102: '(n,gamma)',
103: '(n,p)', 104: '(n,d)', 105: '(n,t)', 106: '(n,3He)',
107: '(n,a)', 108: '(n,2a)', 109: '(n,3a)', 111: '(n,2p)',
112: '(n,pa)', 113: '(n,t2a)', 114: '(n,d2a)', 115: '(n,pd)',
116: '(n,pt)', 117: '(n,da)', 152: '(n,5n)', 153: '(n,6n)',
154: '(n,2nt)', 155: '(n,ta)', 156: '(n,4np)', 157: '(n,3nd)',
158: '(n,nda)', 159: '(n,2npa)', 160: '(n,7n)', 161: '(n,8n)',
162: '(n,5np)', 163: '(n,6np)', 164: '(n,7np)', 165: '(n,4na)',
166: '(n,5na)', 167: '(n,6na)', 168: '(n,7na)', 169: '(n,4nd)',
170: '(n,5nd)', 171: '(n,6nd)', 172: '(n,3nt)', 173: '(n,4nt)',
174: '(n,5nt)', 175: '(n,6nt)', 176: '(n,2n3He)',
177: '(n,3n3He)', 178: '(n,4n3He)', 179: '(n,3n2p)',
180: '(n,3n3a)', 181: '(n,3npa)', 182: '(n,dt)',
183: '(n,npd)', 184: '(n,npt)', 185: '(n,ndt)',
186: '(n,np3He)', 187: '(n,nd3He)', 188: '(n,nt3He)',
189: '(n,nta)', 190: '(n,2n2p)', 191: '(n,p3He)',
192: '(n,d3He)', 193: '(n,3Hea)', 194: '(n,4n2p)',
195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)'}
REACTION_NAME.update({i: '(n,n{})'.format(i-50) for i in range(50, 91)})
REACTION_NAME.update({i: '(n,p{})'.format(i-600) for i in range(600, 649)})
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)})
def _get_products(ev, mt):
"""Generate products from MF=6 in an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation to read from
mt : int
The MT value of the reaction to get products for
Returns
-------
products : list of openmc.data.Product
Products of the reaction
"""
file_obj = StringIO(ev.section[6, mt])
# Read HEAD record
items = get_head_record(file_obj)
reference_frame = {1: 'laboratory', 2: 'center-of-mass',
3: 'light-heavy'}[items[3]]
n_products = items[4]
products = []
for i in range(n_products):
# Get yield for this product
params, yield_ = get_tab1_record(file_obj)
za = params[0]
awr = params[1]
lip = params[2]
law = params[3]
if za == 0:
p = Product('photon')
elif za == 1:
p = Product('neutron')
elif za == 1000:
p = Product('electron')
else:
z = za // 1000
a = za % 1000
p = Product('{}{}'.format(ATOMIC_SYMBOL[z], a))
p.yield_ = yield_
"""
# Set reference frame
if reference_frame == 'laboratory':
p.center_of_mass = False
elif reference_frame == 'center-of-mass':
p.center_of_mass = True
elif reference_frame == 'light-heavy':
p.center_of_mass = (awr <= 4.0)
"""
if law == 0:
# No distribution given
pass
if law == 1:
# Continuum energy-angle distribution
# Peak ahead to determine type of distribution
position = file_obj.tell()
params = get_cont_record(file_obj)
file_obj.seek(position)
lang = params[2]
if lang == 1:
p.distribution = [CorrelatedAngleEnergy.from_endf(file_obj)]
elif lang == 2:
p.distribution = [KalbachMann.from_endf(file_obj)]
elif law == 2:
# Discrete two-body scattering
params, tab2 = get_tab2_record(file_obj)
ne = params[5]
energy = np.zeros(ne)
mu = []
for i in range(ne):
items, values = get_list_record(file_obj)
energy[i] = items[1]
lang = items[2]
if lang == 0:
mu.append(Legendre(values))
elif lang == 12:
mu.append(Tabular(values[::2], values[1::2]))
elif lang == 14:
mu.append(Tabular(values[::2], values[1::2],
'log-linear'))
angle_dist = AngleDistribution(energy, mu)
dist = UncorrelatedAngleEnergy(angle_dist)
p.distribution = [dist]
# TODO: Add level-inelastic info?
elif law == 3:
# Isotropic discrete emission
p.distribution = [UncorrelatedAngleEnergy()]
# TODO: Add level-inelastic info?
elif law == 4:
# Discrete two-body recoil
pass
elif law == 5:
# Charged particle elastic scattering
pass
elif law == 6:
# N-body phase-space distribution
p.distribution = [NBodyPhaseSpace.from_endf(file_obj)]
elif law == 7:
# Laboratory energy-angle distribution
p.distribution = [LaboratoryAngleEnergy.from_endf(file_obj)]
products.append(p)
return products
def _get_fission_products_ace(ace):
"""Generate fission products from an ACE table
Parameters
@ -148,7 +317,116 @@ def _get_fission_products(ace):
return products, derived_products
def _get_photon_products(ace, rx):
def _get_fission_products_endf(ev):
"""Generate fission products from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
Returns
-------
products : list of openmc.data.Product
Prompt and delayed fission neutrons
derived_products : list of openmc.data.Product
"Total" fission neutron
"""
products = []
derived_products = []
if (1, 456) in ev.section:
prompt_neutron = Product('neutron')
prompt_neutron.emission_mode = 'prompt'
# Prompt nu values
file_obj = StringIO(ev.section[1, 456])
lnu = get_head_record(file_obj)[3]
if lnu == 1:
# Polynomial representation
items, coefficients = get_list_record(file_obj)
prompt_neutron.yield_ = Polynomial(coefficients)
elif lnu == 2:
# Tabulated representation
params, prompt_neutron.yield_ = get_tab1_record(file_obj)
products.append(prompt_neutron)
if (1, 452) in ev.section:
total_neutron = Product('neutron')
total_neutron.emission_mode = 'total'
# Total nu values
file_obj = StringIO(ev.section[1, 452])
lnu = get_head_record(file_obj)[3]
if lnu == 1:
# Polynomial representation
items, coefficients = get_list_record(file_obj)
total_neutron.yield_ = Polynomial(coefficients)
elif lnu == 2:
# Tabulated representation
params, total_neutron.yield_ = get_tab1_record(file_obj)
if (1, 456) in ev.section:
derived_products.append(total_neutron)
else:
products.append(total_neutron)
if (1, 455) in ev.section:
file_obj = StringIO(ev.section[1, 455])
# Determine representation of delayed nu data
items = get_head_record(file_obj)
ldg = items[2]
lnu = items[3]
if ldg == 0:
# Delayed-group constants energy independent
items, decay_constants = get_list_record(file_obj)
for constant in decay_constants:
delayed_neutron = Product('neutron')
delayed_neutron.emission_mode = 'delayed'
delayed_neutron.decay_rate = constant
products.append(delayed_neutron)
elif ldg == 1:
# Delayed-group constants energy dependent
raise NotImplementedError('Delayed neutron with energy-dependent '
'group constants.')
# In MF=1, MT=455, the delayed-group abundances are actually not
# specified if the group constants are energy-independent. In this case,
# the abundances must be inferred from MF=5, MT=455 where multiple
# energy distributions are given.
if lnu == 1:
# Nu represented as polynomial
items, coefficients = get_list_record(file_obj)
yield_ = Polynomial(coefficients)
for neutron in products[-6:]:
neutron.yield_ = deepcopy(yield_)
elif lnu == 2:
# Nu represented by tabulation
params, yield_ = get_tab1_record(file_obj)
for neutron in products[-6:]:
neutron.yield_ = deepcopy(yield_)
if (5, 455) in ev.section:
file_obj = StringIO(ev.section[5, 455])
items = get_head_record(file_obj)
nk = items[4]
assert nk == len(decay_constants)
for i in range(nk):
params, applicability = get_tab1_record(file_obj)
dist = UncorrelatedAngleEnergy()
dist.energy = EnergyDistribution.from_endf(file_obj, params)
delayed_neutron = products[-6 + i]
delayed_neutron.yield_.y *= applicability.y[0]
delayed_neutron.distribution.append(dist)
return products, derived_products
def _get_photon_products_ace(ace, rx):
"""Generate photon products from an ACE table
Parameters
@ -253,6 +531,117 @@ def _get_photon_products(ace, rx):
return photons
def _get_photon_products_endf(ev, rx):
"""Generate photon products from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation to read from
rx : openmc.data.Reaction
Reaction that generates photons
Returns
-------
photons : list of openmc.Products
Photons produced from reaction with given MT
"""
products = []
if (12, rx.mt) in ev.section:
file_obj = StringIO(ev.section[12, rx.mt])
items = get_head_record(file_obj)
option = items[2]
if option == 1:
# Multiplicities given
n_discrete_photon = items[4]
if n_discrete_photon > 1:
items, total_yield = get_tab1_record(file_obj)
for k in range(n_discrete_photon):
photon = Product('photon')
# Get photon yield
items, photon.yield_ = get_tab1_record(file_obj)
# Get photon energy distribution
law = items[3]
dist = UncorrelatedAngleEnergy()
if law == 1:
# TODO: Get file 15 distribution
pass
elif law == 2:
energy = items[1]
primary_flag = items[2]
dist.energy = DiscretePhoton(primary_flag, energy,
ev.target['mass'])
photon.distribution.append(dist)
products.append(photon)
elif option == 2:
# Transition probability arrays given
ppyield = {}
ppyield['type'] = 'transition'
ppyield['transition'] = transition = {}
# Determine whether simple (LG=1) or complex (LG=2) transitions
lg = items[3]
# Get transition data
items, values = get_list_record(file_obj)
transition['energy_start'] = items[0]
transition['energies'] = np.array(values[::lg + 1])
transition['direct_probability'] = np.array(values[1::lg + 1])
if lg == 2:
# Complex case
transition['conditional_probability'] = np.array(
values[2::lg + 1])
elif (13, rx.mt) in ev.section:
file_obj = StringIO(ev.section[13, rx.mt])
# Determine option
items = get_head_record(file_obj)
n_discrete_photon = items[4]
if n_discrete_photon > 1:
items, total_xs = get_tab1_record(file_obj)
for k in range(n_discrete_photon):
photon = Product('photon')
items, xs = get_tab1_record(file_obj)
# Re-interpolation photon production cross section and neutron cross
# section to union energy grid
energy = np.union1d(xs.x, rx.xs['0K'].x)
photon_prod_xs = xs(energy)
neutron_xs = rx.xs['0K'](energy)
idx = np.where(neutron_xs > 0)
# Calculate yield as ratio
yield_ = np.zeros_like(energy)
yield_[idx] = photon_prod_xs[idx] / neutron_xs[idx]
photon.yield_ = Tabulated1D(energy, yield_)
# Get photon energy distribution
law = items[3]
dist = UncorrelatedAngleEnergy()
if law == 1:
# TODO: Get file 15 distribution
pass
elif law == 2:
energy = items[1]
primary_flag = items[2]
dist.energy = DiscretePhoton(primary_flag, energy,
ev.target['mass'])
photon.distribution.append(dist)
products.append(photon)
return products
class Reaction(EqualityMixin):
"""A nuclear reaction
@ -350,7 +739,7 @@ class Reaction(EqualityMixin):
cv.check_type('reaction cross section dictionary', xs, MutableMapping)
for key, value in xs.items():
cv.check_type('reaction cross section temperature', key, basestring)
cv.check_type('reaction cross section', value, Function1D)
cv.check_type('reaction cross section', value, Callable)
self._xs = xs
def to_hdf5(self, group):
@ -397,7 +786,7 @@ class Reaction(EqualityMixin):
Returns
-------
openmc.data.ace.Reaction
openmc.data.Reaction
Reaction data
"""
@ -500,7 +889,7 @@ class Reaction(EqualityMixin):
rx.products.append(neutron)
else:
assert mt in (18, 19, 20, 21, 38)
rx.products, rx.derived_products = _get_fission_products(ace)
rx.products, rx.derived_products = _get_fission_products_ace(ace)
for p in rx.products:
if p.emission_mode in ('prompt', 'total'):
@ -568,6 +957,92 @@ class Reaction(EqualityMixin):
# ======================================================================
# PHOTON PRODUCTION
rx.products += _get_photon_products(ace, rx)
rx.products += _get_photon_products_ace(ace, rx)
return rx
@classmethod
def from_endf(cls, ev, mt):
"""Generate a reaction from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
mt : int
The MT value of the reaction to get angular distributions for
Returns
-------
rx : openmc.data.Reaction
Reaction data
"""
rx = Reaction(mt)
# Integrated cross section
if (3, mt) in ev.section:
file_obj = StringIO(ev.section[3, mt])
get_head_record(file_obj)
params, rx.xs['0K'] = get_tab1_record(file_obj)
rx.q_value = params[1]
# Get fission product yields (nu) as well as delayed neutron energy
# distributions
if mt in (18, 19, 20, 21, 38):
rx.products, rx.derived_products = _get_fission_products_endf(ev)
if (6, mt) in ev.section:
# Product angle-energy distribution
rx.products = _get_products(ev, mt)
elif (4, mt) in ev.section or (5, mt) in ev.section:
# Uncorrelated angle-energy distribution
neutron = Product('neutron')
# Note that the energy distribution for MT=455 is read in
# _get_fission_products_endf rather than here
if (5, mt) in ev.section:
file_obj = StringIO(ev.section[5, mt])
items = get_head_record(file_obj)
nk = items[4]
for i in range(nk):
params, applicability = get_tab1_record(file_obj)
dist = UncorrelatedAngleEnergy()
dist.energy = EnergyDistribution.from_endf(file_obj, params)
neutron.applicability.append(applicability)
neutron.distribution.append(dist)
elif mt == 2:
# Elastic scattering -- no energy distribution is given since it
# can be calulcated analytically
dist = UncorrelatedAngleEnergy()
neutron.distribution.append(dist)
elif mt >= 51 and mt < 91:
# Level inelastic scattering -- no energy distribution is given
# since it can be calculated analytically. Here we determine the
# necessary parameters to create a LevelInelastic object
dist = UncorrelatedAngleEnergy()
A = ev.target['mass']
threshold = (A + 1.)/A*abs(rx.q_value)
mass_ratio = (A/(A + 1.))**2
dist.energy = LevelInelastic(threshold, mass_ratio)
neutron.distribution.append(dist)
if (4, mt) in ev.section:
for dist in neutron.distribution:
dist.angle = AngleDistribution.from_endf(ev, mt)
if mt in (18, 19, 20, 21, 38) and (5, mt) in ev.section:
# For fission reactions,
rx.products[0].applicability = neutron.applicability
rx.products[0].distribution = neutron.distribution
else:
rx.products.append(neutron)
if (12, mt) in ev.section or (13, mt) in ev.section:
rx.products += _get_photon_products_endf(ev, rx)
return rx

515
openmc/data/reconstruct.pyx Normal file
View file

@ -0,0 +1,515 @@
from libc.stdlib cimport malloc, calloc, free
from libc.math cimport cos, sin, sqrt, atan, M_PI
cimport numpy as np
import numpy as np
from numpy.linalg import inv
cimport cython
cdef extern from "complex.h":
double cabs(double complex)
double complex conj(double complex)
double creal(complex double)
double cimag(complex double)
double complex cexp(double complex)
# Physical constants are from CODATA 2014
cdef double NEUTRON_MASS_ENERGY = 939.5654133e6 # MeV/c^2
cdef double HBAR_C = 197.3269788e5 # MeV-b^0.5
@cython.cdivision(True)
def wave_number(double A, double E):
r"""Neutron wave number in center-of-mass system.
ENDF-102 defines the neutron wave number in the center-of-mass system in
Equation D.10 as
.. math::
k = \frac{2m_n}{\hbar} \frac{A}{A + 1} \sqrt{|E|}
Parameters
----------
A : double
Ratio of target mass to neutron mass
E : double
Energy in eV
Returns
-------
double
Neutron wave number in b^-0.5
"""
return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C
@cython.cdivision(True)
cdef double _wave_number(double A, double E):
return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C
@cython.cdivision(True)
cdef double phaseshift(int l, double rho):
"""Calculate hardsphere phase shift as given in ENDF-102, Equation D.13
Parameters
----------
l : int
Angular momentum quantum number
rho : float
Product of the wave number and the channel radius
Returns
-------
double
Hardsphere phase shift
"""
if l == 0:
return rho
elif l == 1:
return rho - atan(rho)
elif l == 2:
return rho - atan(3*rho/(3 - rho**2))
elif l == 3:
return rho - atan((15*rho - rho**3)/(15 - 6*rho**2))
elif l == 4:
return rho - atan((105*rho - 10*rho**3)/(105 - 45*rho**2 + rho**4))
@cython.cdivision(True)
def penetration_shift(int l, double rho):
r"""Calculate shift and penetration factors as given in ENDF-102, Equations D.11
and D.12.
Parameters
----------
l : int
Angular momentum quantum number
rho : float
Product of the wave number and the channel radius
Returns
-------
double
Penetration factor for given :math:`l`
double
Shift factor for given :math:`l`
"""
cdef double den
if l == 0:
return rho, 0.
elif l == 1:
den = 1 + rho**2
return rho**3/den, -1/den
elif l == 2:
den = 9 + 3*rho**2 + rho**4
return rho**5/den, -(18 + 3*rho**2)/den
elif l == 3:
den = 225 + 45*rho**2 + 6*rho**4 + rho**6
return rho**7/den, -(675 + 90*rho**2 + 6*rho**4)/den
elif l == 4:
den = 11025 + 1575*rho**2 + 135*rho**4 + 10*rho**6 + rho**8
return rho**9/den, -(44100 + 4725*rho**2 + 270*rho**4 + 10*rho**6)/den
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_mlbw(mlbw, double E):
"""Evaluate cross section using MLBW data.
Parameters
----------
mlbw : openmc.data.MultiLevelBreitWigner
Multi-level Breit-Wigner resonance parameters
E : double
Energy in eV at which to evaluate the cross section
Returns
-------
elastic : double
Elastic scattering cross section in barns
capture : double
Radiative capture cross section in barns
fission : double
Fission cross section in barns
"""
cdef int i, nJ, ij, l, n_res, i_res
cdef double elastic, capture, fission
cdef double A, k, rho, rhohat, I
cdef double P, S, phi, cos2phi, sin2phi
cdef double Ex, Q, rhoc, rhochat, P_c, S_c
cdef double jmin, jmax, j, Dl
cdef double E_r, gt, gn, gg, gf, gx, P_r, S_r, P_rx
cdef double gnE, gtE, Eprime, x, f
cdef double *g
cdef double (*s)[2]
cdef double [:,:] params
I = mlbw.target_spin
A = mlbw.atomic_weight_ratio
k = _wave_number(A, E)
elastic = 0.
capture = 0.
fission = 0.
for i, l in enumerate(mlbw._l_values):
params = mlbw._parameter_matrix[l]
rho = k*mlbw.channel_radius[l](E)
rhohat = k*mlbw.scattering_radius[l](E)
P, S = penetration_shift(l, rho)
phi = phaseshift(l, rhohat)
cos2phi = cos(2*phi)
sin2phi = sin(2*phi)
# Determine shift and penetration at modified energy
if mlbw._competitive[i]:
Ex = E + mlbw.q_value[l]*(A + 1)/A
rhoc = mlbw.channel_radius[l](Ex)
rhochat = mlbw.scattering_radius[l](Ex)
P_c, S_c = penetration_shift(l, rhoc)
if Ex < 0:
P_c = 0
# Determine range of total angular momentum values based on equation
# 41 in LA-UR-12-27079
jmin = abs(abs(I - l) - 0.5)
jmax = I + l + 0.5
nJ = int(jmax - jmin + 1)
# Determine Dl factor using Equation 43 in LA-UR-12-27079
Dl = 2*l + 1
g = <double *> malloc(nJ*sizeof(double))
for ij in range(nJ):
j = jmin + ij
g[ij] = (2*j + 1)/(4*I + 2)
Dl -= g[ij]
s = <double (*)[2]> calloc(2*nJ, sizeof(double))
for i_res in range(params.shape[0]):
# Copy resonance parameters
E_r = params[i_res, 0]
j = params[i_res, 2]
ij = int(j - jmin)
gt = params[i_res, 3]
gn = params[i_res, 4]
gg = params[i_res, 5]
gf = params[i_res, 6]
gx = params[i_res, 7]
P_r = params[i_res, 8]
S_r = params[i_res, 9]
P_rx = params[i_res, 10]
# Calculate neutron and total width at energy E
gnE = P*gn/P_r # ENDF-102, Equation D.7
gtE = gnE + gg + gf
if gx > 0:
gtE += gx*P_c/P_rx
Eprime = E_r + (S_r - S)/(2*P_r)*gn # ENDF-102, Equation D.9
x = 2*(E - Eprime)/gtE # LA-UR-12-27079, Equation 26
f = 2*gnE/(gtE*(1 + x*x)) # Common factor in Equation 40
s[ij][0] += f # First sum in Equation 40
s[ij][1] += f*x # Second sum in Equation 40
capture += f*g[ij]*gg/gtE
if gf > 0:
fission += f*g[ij]*gf/gtE
for ij in range(nJ):
# Add all but last term of LA-UR-12-27079, Equation 40
elastic += g[ij]*((1 - cos2phi - s[ij][0])**2 +
(sin2phi + s[ij][1])**2)
# Add final term with Dl from Equation 40
elastic += 2*Dl*(1 - cos2phi)
# Free memory
free(g)
free(s)
capture *= 2*M_PI/(k*k)
fission *= 2*M_PI/(k*k)
elastic *= M_PI/(k*k)
return (elastic, capture, fission)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_slbw(slbw, double E):
"""Evaluate cross section using SLBW data.
Parameters
----------
slbw : openmc.data.SingleLevelBreitWigner
Single-level Breit-Wigner resonance parameters
E : double
Energy in eV at which to evaluate the cross section
Returns
-------
elastic : double
Elastic scattering cross section in barns
capture : double
Radiative capture cross section in barns
fission : double
Fission cross section in barns
"""
cdef int i, l, i_res
cdef double elastic, capture, fission
cdef double A, k, rho, rhohat, I
cdef double P, S, phi, cos2phi, sin2phi, sinphi2
cdef double Ex, rhoc, rhochat, P_c, S_c
cdef double E_r, J, gt, gn, gg, gf, gx, P_r, S_r, P_rx
cdef double gnE, gtE, Eprime, f
cdef double x, theta, psi, chi
cdef double [:,:] params
I = slbw.target_spin
A = slbw.atomic_weight_ratio
k = _wave_number(A, E)
elastic = 0.
capture = 0.
fission = 0.
for i, l in enumerate(slbw._l_values):
params = slbw._parameter_matrix[l]
rho = k*slbw.channel_radius[l](E)
rhohat = k*slbw.scattering_radius[l](E)
P, S = penetration_shift(l, rho)
phi = phaseshift(l, rhohat)
cos2phi = cos(2*phi)
sin2phi = sin(2*phi)
sinphi2 = sin(phi)**2
# Add potential scattering -- first term in ENDF-102, Equation D.2
elastic += 4*M_PI/(k*k)*(2*l + 1)*sinphi2
# Determine shift and penetration at modified energy
if slbw._competitive[i]:
Ex = E + slbw.q_value[l]*(A + 1)/A
rhoc = slbw.channel_radius[l](Ex)
rhochat = slbw.scattering_radius[l](Ex)
P_c, S_c = penetration_shift(l, rhoc)
if Ex < 0:
P_c = 0
for i_res in range(params.shape[0]):
# Copy resonance parameters
E_r = params[i_res, 0]
J = params[i_res, 2]
gt = params[i_res, 3]
gn = params[i_res, 4]
gg = params[i_res, 5]
gf = params[i_res, 6]
gx = params[i_res, 7]
P_r = params[i_res, 8]
S_r = params[i_res, 9]
P_rx = params[i_res, 10]
# Calculate neutron and total width at energy E
gnE = P*gn/P_r # Equation D.7
gtE = gnE + gg + gf
if gx > 0:
gtE += gx*P_c/P_rx
Eprime = E_r + (S_r - S)/(2*P_r)*gn # Equation D.9
gJ = (2*J + 1)/(4*I + 2) # Mentioned in section D.1.1.4
# Calculate common factor for elastic, capture, and fission
# cross sections
f = M_PI/(k*k)*gJ*gnE/((E - Eprime)**2 + gtE**2/4)
# Add contribution to elastic per Equation D.2
elastic += f*(gnE*cos2phi - 2*(gg + gf)*sinphi2
+ 2*(E - Eprime)*sin2phi)
# Add contribution to capture per Equation D.3
capture += f*gg
# Add contribution to fission per Equation D.6
if gf > 0:
fission += f*gf
return (elastic, capture, fission)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_rm(rm, double E):
"""Evaluate cross section using Reich-Moore data.
Parameters
----------
rm : openmc.data.ReichMoore
Reich-Moore resonance parameters
E : double
Energy in eV at which to evaluate the cross section
Returns
-------
elastic : double
Elastic scattering cross section in barns
capture : double
Radiative capture cross section in barns
fission : double
Fission cross section in barns
"""
cdef int i, l, m, n, i_res
cdef int i_s, num_s, i_J, num_J
cdef double elastic, capture, fission, total
cdef double A, k, rho, rhohat, I
cdef double P, S, phi
cdef double smin, smax, s, Jmin, Jmax, J, j
cdef double E_r, gn, gg, gfa, gfb, P_r
cdef double E_diff, abs_value, gJ
cdef double complex Ubar, U_, denominator_inverse
cdef bint hasfission
cdef np.ndarray[double, ndim=2] one
cdef np.ndarray[double complex, ndim=2] K, Imat, U
cdef double [:,:] params
# Get nuclear spin
I = rm.target_spin
elastic = 0.
fission = 0.
total = 0.
A = rm.atomic_weight_ratio
k = _wave_number(A, E)
one = np.eye(3)
K = np.zeros((3,3), dtype=complex)
for i, l in enumerate(rm._l_values):
# Check for l-dependent scattering radius
rho = k*rm.channel_radius[l](E)
rhohat = k*rm.scattering_radius[l](E)
# Calculate shift and penetrability
P, S = penetration_shift(l, rho)
# Calculate phase shift
phi = phaseshift(l, rhohat)
# Calculate common factor on collision matrix terms (term outside curly
# braces in ENDF-102, Eq. D.27)
Ubar = cexp(-2j*phi)
# The channel spin is the vector sum of the target spin, I, and the
# neutron spin, 1/2, so can take on values of |I - 1/2| < s < I + 1/2
smin = abs(I - 0.5)
smax = I + 0.5
num_s = int(smax - smin + 1)
for i_s in range(num_s):
s = i_s + smin
# Total angular momentum is the vector sum of l and s and can assume
# values between |l - s| < J < l + s
Jmin = abs(l - s)
Jmax = l + s
num_J = int(Jmax - Jmin + 1)
for i_J in range(num_J):
J = i_J + Jmin
# Initialize K matrix
for m in range(3):
for n in range(3):
K[m,n] = 0.0
hasfission = False
if (l, J) in rm._parameter_matrix:
params = rm._parameter_matrix[l, J]
for i_res in range(params.shape[0]):
# Sometimes, the same (l, J) quantum numbers can occur
# for different values of the channel spin, s. In this
# case, the sign of the channel spin indicates which
# spin is to be used. If the spin is negative assume
# this resonance comes from the I - 1/2 channel and vice
# versa.
j = params[i_res, 2]
if l > 0:
if (j < 0 and s != smin) or (j > 0 and s != smax):
continue
# Copy resonance parameters
E_r = params[i_res, 0]
gn = params[i_res, 3]
gg = params[i_res, 4]
gfa = params[i_res, 5]
gfb = params[i_res, 6]
P_r = params[i_res, 7]
# Calculate neutron width at energy E
gn = sqrt(P*gn/P_r)
# Calculate inverse of denominator of K matrix terms
E_diff = E_r - E
abs_value = E_diff*E_diff + 0.25*gg*gg
denominator_inverse = (E_diff + 0.5j*gg)/abs_value
# Upper triangular portion of K matrix -- see ENDF-102,
# Equation D.28
K[0,0] = K[0,0] + gn*gn*denominator_inverse
if gfa != 0.0 or gfb != 0.0:
# Negate fission widths if necessary
gfa = (-1 if gfa < 0 else 1)*sqrt(abs(gfa))
gfb = (-1 if gfb < 0 else 1)*sqrt(abs(gfb))
K[0,1] = K[0,1] + gn*gfa*denominator_inverse
K[0,2] = K[0,2] + gn*gfb*denominator_inverse
K[1,1] = K[1,1] + gfa*gfa*denominator_inverse
K[1,2] = K[1,2] + gfa*gfb*denominator_inverse
K[2,2] = K[2,2] + gfb*gfb*denominator_inverse
hasfission = True
# multiply by factor of i/2
for m in range(3):
for n in range(3):
K[m,n] = K[m,n] * 0.5j
# Get collision matrix
gJ = (2*J + 1)/(4*I + 2)
if hasfission:
# Copy upper triangular portion of K to lower triangular
K[1,0] = K[0,1]
K[2,0] = K[0,2]
K[2,1] = K[1,2]
Imat = inv(one - K)
U = Ubar*(2*Imat - one) # ENDF-102, Eq. D.27
elastic += gJ*cabs(1 - U[0,0])**2 # ENDF-102, Eq. D.24
total += 2*gJ*(1 - creal(U[0,0])) # ENDF-102, Eq. D.23
# Calculate fission from ENDF-102, Eq. D.26
fission += 4*gJ*(cabs(Imat[1,0])**2 + cabs(Imat[2,0])**2)
else:
U_ = Ubar*(2*conj(1 - K[0,0])/cabs(1 - K[0,0])**2 - 1)
elastic += gJ*cabs(1 - U_)**2 # ENDF-102, Eq. D.24
total += 2*gJ*(1 - creal(U_)) # ENDF-102, Eq. D.23
# Calculate capture as difference of other cross sections as per ENDF-102,
# Equation D.25
capture = total - elastic - fission
elastic *= M_PI/(k*k)
capture *= M_PI/(k*k)
fission *= M_PI/(k*k)
return (elastic, capture, fission)

1059
openmc/data/resonance.py Normal file
View file

File diff suppressed because it is too large Load diff

12
setup.py
View file

@ -8,6 +8,12 @@ except ImportError:
from distutils.core import setup
have_setuptools = False
try:
from Cython.Build import cythonize
have_cython = True
except ImportError:
have_cython = False
kwargs = {'name': 'openmc',
'version': '0.8.0',
'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.model',
@ -48,4 +54,10 @@ if have_setuptools:
},
})
# If Cython is present, add resonance reconstruction capability
if have_cython:
kwargs.update({
'ext_modules': cythonize('openmc/data/reconstruct.pyx')
})
setup(**kwargs)