From cbcd3a7f573c1c5ead8a7d1ec5457d34fa034f87 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Feb 2016 11:01:12 -0600 Subject: [PATCH] Extend openmc.data to be able to import ENDF data and perform resonance reconstruction. --- .gitignore | 5 + docs/source/pythonapi/index.rst | 73 +- openmc/data/__init__.py | 2 + openmc/data/angle_distribution.py | 99 ++- openmc/data/angle_energy.py | 3 +- openmc/data/correlated.py | 53 +- openmc/data/data.py | 60 +- openmc/data/endf.py | 409 +++++++++++ openmc/data/endf_utils.py | 44 -- openmc/data/energy_distribution.py | 189 +++++ openmc/data/fission_energy.py | 68 +- openmc/data/function.py | 78 ++ openmc/data/kalbach_mann.py | 52 ++ openmc/data/laboratory.py | 139 ++++ openmc/data/nbody.py | 23 + openmc/data/neutron.py | 100 ++- openmc/data/product.py | 3 +- openmc/data/reaction.py | 493 ++++++++++++- openmc/data/reconstruct.pyx | 515 ++++++++++++++ openmc/data/resonance.py | 1059 ++++++++++++++++++++++++++++ setup.py | 12 + 21 files changed, 3308 insertions(+), 171 deletions(-) create mode 100644 openmc/data/endf.py delete mode 100644 openmc/data/endf_utils.py create mode 100644 openmc/data/laboratory.py create mode 100644 openmc/data/reconstruct.pyx create mode 100644 openmc/data/resonance.py diff --git a/.gitignore b/.gitignore index 2ca1c7d14..b72a18c84 100644 --- a/.gitignore +++ b/.gitignore @@ -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 diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index d8eb5200a..93507442e 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -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/ diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index e60878d60..3f42ef90e 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -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 * diff --git a/openmc/data/angle_distribution.py b/openmc/data/angle_distribution.py index 80b6cdbdb..0f0105c2a 100644 --- a/openmc/data/angle_distribution.py +++ b/openmc/data/angle_distribution.py @@ -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) diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 4a5e391e3..9cee89aa5 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -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 ---------- diff --git a/openmc/data/correlated.py b/openmc/data/correlated.py index d6ce8ab47..dac1968ad 100644 --- a/openmc/data/correlated.py +++ b/openmc/data/correlated.py @@ -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) diff --git a/openmc/data/data.py b/openmc/data/data.py index 574ae541a..46fcdf36f 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -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]) diff --git a/openmc/data/endf.py b/openmc/data/endf.py new file mode 100644 index 000000000..ba398598d --- /dev/null +++ b/openmc/data/endf.py @@ -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 diff --git a/openmc/data/endf_utils.py b/openmc/data/endf_utils.py deleted file mode 100644 index 1a77c60a5..000000000 --- a/openmc/data/endf_utils.py +++ /dev/null @@ -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} diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py index a16090635..b42887b2f 100644 --- a/openmc/data/energy_distribution.py +++ b/openmc/data/energy_distribution.py @@ -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 diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 6352915cc..c7f71b657 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -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 diff --git a/openmc/data/function.py b/openmc/data/function.py index 94a2f0911..13df12b47 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -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 diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index 5aca17ba9..e84e6a77a 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -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) diff --git a/openmc/data/laboratory.py b/openmc/data/laboratory.py new file mode 100644 index 000000000..be449b79a --- /dev/null +++ b/openmc/data/laboratory.py @@ -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) diff --git a/openmc/data/nbody.py b/openmc/data/nbody.py index a34380dc9..6616b1858 100644 --- a/openmc/data/nbody.py +++ b/openmc/data/nbody.py @@ -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) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 2431ba74a..0c74aeaa1 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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 diff --git a/openmc/data/product.py b/openmc/data/product.py index 753888f6d..d9c3506e5 100644 --- a/openmc/data/product.py +++ b/openmc/data/product.py @@ -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 diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index eb33e1395..cef635135 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -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 diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx new file mode 100644 index 000000000..4790cc345 --- /dev/null +++ b/openmc/data/reconstruct.pyx @@ -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 = malloc(nJ*sizeof(double)) + for ij in range(nJ): + j = jmin + ij + g[ij] = (2*j + 1)/(4*I + 2) + Dl -= g[ij] + + s = 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) diff --git a/openmc/data/resonance.py b/openmc/data/resonance.py new file mode 100644 index 000000000..206a56df4 --- /dev/null +++ b/openmc/data/resonance.py @@ -0,0 +1,1059 @@ +from collections import defaultdict, MutableSequence, Iterable +import io + +import numpy as np +from numpy.polynomial import Polynomial +import pandas as pd + +from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record +try: + from .reconstruct import wave_number, penetration_shift, reconstruct_mlbw, \ + reconstruct_slbw, reconstruct_rm + _reconstruct = True +except ImportError: + _reconstruct = False +import openmc.checkvalue as cv + +class Resonances(object): + """Resolved and unresolved resonance data + + Parameters + ---------- + ranges : list of openmc.data.ResonanceRange + Distinct energy ranges for resonance data + + Attributes + ---------- + ranges : list of openmc.data.ResonanceRange + Distinct energy ranges for resonance data + resolved : openmc.data.ResonanceRange or None + Resolved resonance range + unresolved : openmc.data.Unresolved or None + Unresolved resonance range + + """ + + def __init__(self, ranges): + self.ranges = ranges + + @property + def ranges(self): + return self._ranges + + @property + def resolved(self): + resolved_ranges = [r for r in self.ranges + if not isinstance(r, Unresolved)] + if len(resolved_ranges) > 1: + raise ValueError('More than one resolved range present') + elif len(resolved_ranges) == 0: + return None + else: + return resolved_ranges[0] + + @property + def unresolved(self): + for r in self.ranges: + if isinstance(r, Unresolved): + return r + else: + return None + + @ranges.setter + def ranges(self, ranges): + cv.check_type('resonance ranges', ranges, MutableSequence) + self._ranges = cv.CheckedList(ResonanceRange, 'resonance ranges', + ranges) + + @classmethod + def from_endf(cls, ev): + """Generate resonance data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + + Returns + ------- + openmc.data.Resonances + Resonance data + + """ + file_obj = io.StringIO(ev.section[2, 151]) + + # Determine whether discrete or continuous representation + items = get_head_record(file_obj) + n_isotope = items[4] # Number of isotopes + + ranges = [] + for iso in range(n_isotope): + items = get_cont_record(file_obj) + abundance = items[1] + fission_widths = (items[3] == 1) # fission widths are given? + n_ranges = items[4] # number of resonance energy ranges + + for j in range(n_ranges): + items = get_cont_record(file_obj) + resonance_flag = items[2] # flag for resolved (1)/unresolved (2) + formalism = items[3] # resonance formalism + + if resonance_flag in (0, 1): + # resolved resonance region + erange = _FORMALISMS[formalism].from_endf(ev, file_obj, items) + + elif resonance_flag == 2: + # unresolved resonance region + erange = Unresolved.from_endf(file_obj, items, fission_widths) + + #erange.material = self + ranges.append(erange) + + return cls(ranges) + + +class ResonanceRange(object): + """Resolved resonance range + + Parameters + ---------- + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + channel : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + scattering : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + + Attributes + ---------- + channel_radius : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + energy_max : float + Maximum energy of the resolved resonance range in eV + energy_min : float + Minimum energy of the resolved resonance range in eV + scattering_radius : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energ + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + + """ + def __init__(self, target_spin, energy_min, energy_max, channel, scattering): + self.target_spin = target_spin + self.energy_min = energy_min + self.energy_max = energy_max + self.channel_radius = channel + self.scattering_radius = scattering + + self._prepared = False + self._parameter_matrix = {} + + @classmethod + def from_endf(cls, ev, file_obj, items): + """Create resonance range from an ENDF evaluation. + + This factory method is only used when LRU=0, indicating that only a + scattering radius appears in MF=2, MT=151. All subclasses of + ResonanceRange override this method with their own. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + + Returns + ------- + openmc.data.ResonanceRange + Resonance range data + + """ + energy_min, energy_max = items[0:2] + + # For scattering radius-only, NRO must be zero + assert items[4] == 0 + + # Get energy-independent scattering radius + items = get_cont_record(file_obj) + target_spin = items[0] + ap = Polynomial((items[1],)) + + # Calculate channel radius from ENDF-102 equation D.14 + a = Polynomial((0.123 * (1.008665*ev.target['mass'])**(1./3.) + 0.08,)) + + return cls(target_spin, energy_min, energy_max, {0: a}, {0: ap}) + + def reconstruct(self, energies): + """Evaluate cross section at specified energies. + + Parameters + ---------- + energies : float or Iterable of float + Energies at which the cross section should be evaluated + + Returns + ------- + 3-tuple of float or numpy.ndarray + Elastic, capture, and fission cross sections at the specified + energies + + """ + if not _reconstruct: + raise RuntimeError("Resonance reconstruction not available.") + + # Pre-calculate penetrations and shifts for resonances + if not self._prepared: + self._prepare_resonances() + + if isinstance(energies, Iterable): + elastic = np.zeros_like(energies) + capture = np.zeros_like(energies) + fission = np.zeros_like(energies) + + for i, E in enumerate(energies): + xse, xsg, xsf = self._reconstruct(self, E) + elastic[i] = xse + capture[i] = xsg + fission[i] = xsf + else: + elastic, capture, fission = self._reconstruct(self, energies) + + return {2: elastic, 102: capture, 18: fission} + + +class MultiLevelBreitWigner(ResonanceRange): + """Multi-level Breit-Wigner resolved resonance formalism data. + + Multi-level Breit-Wigner resolved resonance data is identified by LRF=2 in + the ENDF-6 format. + + Parameters + ---------- + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + channel : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + scattering : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + + Attributes + ---------- + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide given as a function of + l-value. Note that this may be different than the value for the + evaluation as a whole. + channel_radius : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + energy_max : float + Maximum energy of the resolved resonance range in eV + energy_min : float + Minimum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Energies, spins, and resonances widths for each resonance + q_value : dict + Q-value to be added to incident particle's center-of-mass energy to + determine the channel energy for use in the penetrability factor. The + keys of the dictionary are l-values. + scattering_radius : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + + """ + + def __init__(self, target_spin, energy_min, energy_max, channel, scattering): + super(MultiLevelBreitWigner, self).__init__( + target_spin, energy_min, energy_max, channel, scattering) + self.parameters = None + self.q_value = {} + self.atomic_weight_ratio = None + + # Set resonance reconstruction function + if _reconstruct: + self._reconstruct = reconstruct_mlbw + else: + self._reconstruct = None + + @classmethod + def from_endf(cls, ev, file_obj, items): + """Create MLBW data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + + Returns + ------- + openmc.data.MultiLevelBreitWigner + Multi-level Breit-Wigner resonance parameters + + """ + + # Read energy-dependent scattering radius if present + energy_min, energy_max = items[0:2] + nro, naps = items[4:6] + if nro != 0: + params, ape = get_tab1_record(file_obj) + + # Other scatter radius parameters + items = get_cont_record(file_obj) + target_spin = items[0] + ap = Polynomial((items[1],)) # energy-independent scattering-radius + NLS = items[4] # number of l-values + + # Read resonance widths, J values, etc + channel_radius = {} + scattering_radius = {} + q_value = {} + records = [] + for l in range(NLS): + items, values = get_list_record(file_obj) + l_value = items[2] + awri = items[0] + q_value[l_value] = items[1] + competitive = items[3] + + # Calculate channel radius from ENDF-102 equation D.14 + a = Polynomial((0.123 * (1.008665*awri)**(1./3.) + 0.08,)) + + # Construct scattering and channel radius + if nro == 0: + scattering_radius[l_value] = ap + if naps == 0: + channel_radius[l_value] = a + elif naps == 1: + channel_radius[l_value] = ap + elif nro == 1: + scattering_radius[l_value] = ape + if naps == 0: + channel_radius[l_value] = a + elif naps == 1: + channel_radius[l_value] = ape + elif naps == 2: + channel_radius[l_value] = ap + + energy = values[0::6] + spin = values[1::6] + gt = np.asarray(values[2::6]) + gn = np.asarray(values[3::6]) + gg = np.asarray(values[4::6]) + gf = np.asarray(values[5::6]) + if competitive > 0: + gx = gt - (gn + gg + gf) + else: + gx = np.zeros_like(gt) + + for i, E in enumerate(energy): + records.append([energy[i], l_value, spin[i], gt[i], gn[i], + gg[i], gf[i], gx[i]]) + + columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', + 'captureWidth', 'fissionWidth', 'competitiveWidth'] + parameters = pd.DataFrame.from_records(records, columns=columns) + + # Create instance of class + mlbw = cls(target_spin, energy_min, energy_max, + channel_radius, scattering_radius) + mlbw.q_value = q_value + mlbw.atomic_weight_ratio = awri + mlbw.parameters = parameters + + return mlbw + + def _prepare_resonances(self): + df = self.parameters.copy() + + # Penetration and shift factors + p = np.zeros(len(df)) + s = np.zeros(len(df)) + + # Penetration and shift factors for competitive reaction + px = np.zeros(len(df)) + sx = np.zeros(len(df)) + + l_values = [] + competitive = [] + + A = self.atomic_weight_ratio + for i, E, l, J, gt, gn, gg, gf, gx in df.itertuples(): + if l not in l_values: + l_values.append(l) + competitive.append(gx > 0) + + # Determine penetration and shift corresponding to resonance energy + k = wave_number(A, E) + rho = k*self.channel_radius[l](E) + rhohat = k*self.scattering_radius[l](E) + p[i], s[i] = penetration_shift(l, rho) + + # Determine penetration at modified energy for competitive reaction + if gx > 0: + Ex = E + self.q[l]*(A + 1)/A + rho = k*self.channel_radius[l](Ex) + rhohat = k*self.scattering_radius[l](Ex) + px[i], sx[i] = penetration_shift(l, rho) + else: + px[i] = sx[i] = 0.0 + + df['p'] = p + df['s'] = s + df['px'] = px + df['sx'] = sx + + self._l_values = np.array(l_values) + self._competitive = np.array(competitive) + for l in l_values: + self._parameter_matrix[l] = df[df.L == l].as_matrix() + + self._prepared = True + + +class SingleLevelBreitWigner(MultiLevelBreitWigner): + """Single-level Breit-Wigner resolved resonance formalism data. + + Single-level Breit-Wigner resolved resonance data is is identified by LRF=1 + in the ENDF-6 format. + + Parameters + ---------- + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + channel : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + scattering : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + + Attributes + ---------- + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide given as a function of + l-value. Note that this may be different than the value for the + evaluation as a whole. + channel_radius : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + energy_max : float + Maximum energy of the resolved resonance range in eV + energy_min : float + Minimum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Energies, spins, and resonances widths for each resonance + q_value : dict + Q-value to be added to incident particle's center-of-mass energy to + determine the channel energy for use in the penetrability factor. The + keys of the dictionary are l-values. + scattering_radius : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + + """ + + def __init__(self, target_spin, energy_min, energy_max, channel, scattering): + super(SingleLevelBreitWigner, self).__init__( + target_spin, energy_min, energy_max, channel, scattering) + + # Set resonance reconstruction function + if _reconstruct: + self._reconstruct = reconstruct_slbw + else: + self._reconstruct = None + + +class ReichMoore(ResonanceRange): + """Reich-Moore resolved resonance formalism data. + + Reich-Moore resolved resonance data is identified by LRF=3 in the ENDF-6 + format. + + Parameters + ---------- + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + channel : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + scattering : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + + Attributes + ---------- + angle_distribution : bool + Indicate whether parameters can be used to compute angular distributions + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide given as a function of + l-value. Note that this may be different than the value for the + evaluation as a whole. + channel_radius : dict + Dictionary whose keys are l-values and values are channel radii as a + function of energy + energy_max : float + Maximum energy of the resolved resonance range in eV + energy_min : float + Minimum energy of the resolved resonance range in eV + num_l_convergence : int + Number of l-values which must be used to converge the calculation + scattering_radius : dict + Dictionary whose keys are l-values and values are scattering radii as a + function of energy + parameters : pandas.DataFrame + Energies, spins, and resonances widths for each resonance + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + + """ + + def __init__(self, target_spin, energy_min, energy_max, channel, scattering): + super(ReichMoore, self).__init__( + target_spin, energy_min, energy_max, channel, scattering) + self.parameters = None + self.angle_distribution = False + self.num_l_convergence = 0 + + # Set resonance reconstruction function + if _reconstruct: + self._reconstruct = reconstruct_rm + else: + self._reconstruct = None + + @classmethod + def from_endf(cls, ev, file_obj, items): + """Create Reich-Moore resonance data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + + Returns + ------- + openmc.data.ReichMoore + Reich-Moore resonance parameters + + """ + # Read energy-dependent scattering radius if present + energy_min, energy_max = items[0:2] + nro, naps = items[4:6] + if nro != 0: + params, ape = get_tab1_record(file_obj) + + # Other scatter radius parameters + items = get_cont_record(file_obj) + target_spin = items[0] + ap = Polynomial((items[1],)) + angle_distribution = (items[3] == 1) # Flag for angular distribution + NLS = items[4] # Number of l-values + num_l_convergence = items[5] # Number of l-values for convergence + + # Read resonance widths, J values, etc + channel_radius = {} + scattering_radius = {} + records = [] + for i in range(NLS): + items, values = get_list_record(file_obj) + apl = Polynomial((items[1],)) if items[1] != 0.0 else ap + l_value = items[2] + awri = items[0] + + # Calculate channel radius from ENDF-102 equation D.14 + a = Polynomial((0.123 * (1.008665*awri)**(1./3.) + 0.08,)) + + # Construct scattering and channel radius + if nro == 0: + scattering_radius[l_value] = apl + if naps == 0: + channel_radius[l_value] = a + elif naps == 1: + channel_radius[l_value] = apl + elif nro == 1: + if naps == 0: + channel_radius[l_value] = a + scattering_radius[l_value] = ape + elif naps == 1: + channel_radius[l_value] = scattering_radius[l_value] = ape + elif naps == 2: + channel_radius[l_value] = apl + scattering_radius[l_value] = ape + + energy = values[0::6] + spin = values[1::6] + gn = values[2::6] + gg = values[3::6] + gfa = values[4::6] + gfb = values[5::6] + + for i, E in enumerate(energy): + records.append([energy[i], l_value, spin[i], gn[i], gg[i], + gfa[i], gfb[i]]) + + # Create pandas DataFrame with resonance data + columns = ['energy', 'L', 'J', 'neutronWidth', 'captureWidth', + 'fissionWidthA', 'fissionWidthB'] + parameters = pd.DataFrame.from_records(records, columns=columns) + + # Create instance of ReichMoore + rm = cls(target_spin, energy_min, energy_max, + channel_radius, scattering_radius) + rm.parameters = parameters + rm.angle_distribution = angle_distribution + rm.num_l_convergence = num_l_convergence + rm.atomic_weight_ratio = awri + + return rm + + def _prepare_resonances(self): + df = self.parameters.copy() + + # Penetration and shift factors + p = np.zeros(len(df)) + s = np.zeros(len(df)) + + l_values = [] + lj_values = [] + + A = self.atomic_weight_ratio + for i, E, l, J, gn, gg, gfa, gfb in df.itertuples(): + if l not in l_values: + l_values.append(l) + if (l, abs(J)) not in lj_values: + lj_values.append((l, abs(J))) + + # Determine penetration and shift corresponding to resonance energy + k = wave_number(A, E) + rho = k*self.channel_radius[l](E) + rhohat = k*self.scattering_radius[l](E) + p[i], s[i] = penetration_shift(l, rho) + + df['p'] = p + df['s'] = s + + self._l_values = np.array(l_values) + for (l, J) in lj_values: + self._parameter_matrix[l, J] = df[(df.L == l) & + (abs(df.J) == J)].as_matrix() + + self._prepared = True + + +class RMatrixLimited(ResonanceRange): + """R-matrix limited resolved resonance formalism data. + + R-matrix limited resolved resonance data is identified by LRF=7 in the + ENDF-6 format. + + Parameters + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + particle_pairs : list of dict + List of particle pairs. Each particle pair is represented by a + dictionary that contains the mass, atomic number, spin, and parity of + each particle as well as other characteristics. + spin_groups : list of dict + List of spin groups. Each spin group is characterized by channels, + resonance energies, and resonance widths. + + Attributes + ---------- + reduced_width : bool + Flag indicating whether channel widths in eV or reduced-width amplitudes + in eV^1/2 are given + formalism : int + Flag to specify which formulae for the R-matrix are to be used + particle_pairs : list of dict + List of particle pairs. Each particle pair is represented by a + dictionary that contains the mass, atomic number, spin, and parity of + each particle as well as other characteristics. + spin_groups : list of dict + List of spin groups. Each spin group is characterized by channels, + resonance energies, and resonance widths. + + """ + + def __init__(self, energy_min, energy_max, particle_pairs, spin_groups): + super(RMatrixLimited, self).__init__(0.0, energy_min, energy_max, + None, None) + self.reduced_width = False + self.formalism = 3 + self.particle_pairs = particle_pairs + self.spin_groups = spin_groups + + @classmethod + def from_endf(cls, ev, file_obj, items): + """Read R-Matrix limited resonance data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + + Returns + ------- + openmc.data.RMatrixLimited + R-matrix limited resonance parameters + + """ + energy_min, energy_max = items[0:2] + + items = get_cont_record(file_obj) + reduced_width = (items[2] == 1) # reduced width amplitude? + formalism = items[3] # Specify which formulae are used + n_spin_groups = items[4] # Number of Jpi values (NJS) + + particle_pairs = [] + spin_groups = [] + + items, values = get_list_record(file_obj) + n_pairs = items[5]//2 # Number of particle pairs (NPP) + for i in range(n_pairs): + first = {'mass': values[12*i], + 'z': int(values[12*i + 2]), + 'spin': values[12*i + 4], + 'parity': values[12*i + 10]} + second = {'mass': values[12*i + 1], + 'z': int(values[12*i + 3]), + 'spin': values[12*i + 5], + 'parity': values[12*i + 11]} + + q_value = values[12*i + 6] + penetrability = values[12*i + 7] + shift = values[12*i + 8] + mt = int(values[12*i + 9]) + + particle_pairs.append(ParticlePair( + first, second, q_value, penetrability, shift, mt)) + + # loop over spin groups + for i in range(n_spin_groups): + items, values = get_list_record(file_obj) + J = items[0] + if J == 0.0: + parity = '+' if items[1] == 1.0 else '-' + else: + parity = '+' if J > 0. else '-' + J = abs(J) + kbk = items[2] + kps = items[3] + n_channels = items[5] + channels = [] + for j in range(n_channels): + channel = {} + channel['particle_pair'] = particle_pairs[ + int(values[6*j]) - 1] + channel['l'] = values[6*j + 1] + channel['spin'] = values[6*j + 2] + channel['boundary'] = values[6*j + 3] + channel['effective_radius'] = values[6*j + 4] + channel['true_radius'] = values[6*j + 5] + channels.append(channel) + + # Read resonance energies and widths + items, values = get_list_record(file_obj) + n_resonances = items[3] + records = [] + m = n_channels//6 + 1 + for j in range(n_resonances): + energy = values[6*m*j] + records.append([energy] + [values[6*m*j + k + 1] + for k in range(n_channels)]) + + # Determine column names + columns = ['energy'] + for channel in channels: + mt = channel['particle_pair'].mt + if mt == 2: + columns.append('neutronWidth') + elif mt == 18: + columns.append('fissionWidth') + elif mt == 102: + columns.append('captureWidth') + else: + columns.append('width (MT={})'.format(mt)) + + # Create Pandas dataframe with resonance parameters + parameters = pd.DataFrame.from_records(records, columns=columns) + + # Construct SpinGroup instance and add to list + sg = SpinGroup(J, parity, channels, parameters) + spin_groups.append(sg) + + # Optional extension (Background R-Matrix) + if kbk > 0: + items, values = get_list_record(file_obj) + lbk = items[4] + if lbk == 1: + params, rbr = get_tab1_record(file_obj) + params, rbi = get_tab1_record(file_obj) + + # Optional extension (Tabulated phase shifts) + if kps > 0: + items, values = get_list_record(file_obj) + lps = items[4] + if lps == 1: + params, psr = get_tab1_record(file_obj) + params, psi = get_tab1_record(file_obj) + + rml = cls(energy_min, energy_max, particle_pairs, spin_groups) + rml.reduced_width = reduced_width + rml.formalism = formalism + + return rml + + +class ParticlePair(object): + def __init__(self, first, second, q_value, penetrability, + shift, mt): + self.first = first + self.second = second + self.q_value = q_value + self.penetrability = penetrability + self.shift = shift + self.mt = mt + + +class SpinGroup(object): + """Resonance spin group + + Attributes + ---------- + spin : float + Total angular momentum (nuclear spin) + parity : {'+', '-'} + Even (+) or odd(-) parity + channels : list of openmc.data.Channel + Available channels + parameters : pandas.DataFrame + Energies/widths for each resonance/channel + + """ + + def __init__(self, spin, parity, channels, parameters): + self.spin = spin + self.parity = parity + self.channels = channels + self.parameters = parameters + + def __repr__(self): + return ''.format(self.spin, self.parity) + + +class Unresolved(ResonanceRange): + """Unresolved resonance parameters as identified by LRU=2 in MF=2. + + Parameters + ---------- + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + energy_min : float + Minimum energy of the unresolved resonance range in eV + energy_max : float + Maximum energy of the unresolved resonance range in eV + scatter : openmc.data.Function1D + Scattering radii as a function of energy + + Attributes + ---------- + add_to_background : bool + If True, file 3 contains partial cross sections to be added to the + average unresolved cross sections calculated from parameters. + energies : Iterable of float + Energies at which parameters are tabulated + energy_max : float + Maximum energy of the unresolved resonance range in eV + energy_min : float + Minimum energy of the unresolved resonance range in eV + parameters : list of pandas.DataFrame + Average resonance parameters at each energy + scattering_radius : openmc.data.Function1D + Scattering radii as a function of energy + target_spin : float + Intrinsic spin, :math:`I`, of the target nuclide + + """ + + def __init__(self, target_spin, energy_min, energy_max, scatter): + super(Unresolved, self).__init__( + target_spin, energy_min, energy_max, None, scatter) + self.energies = None + self.parameters = None + self.add_to_background = False + + @classmethod + def from_endf(cls, file_obj, items, fission_widths): + """Read unresolved resonance data from an ENDF evaluation. + + Parameters + ---------- + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + fission_widths : bool + Whether fission widths are given + + Returns + ------- + openmc.data.Unresolved + Unresolved resonance region parameters + + """ + # Read energy-dependent scattering radius if present + energy_min, energy_max = items[0:2] + nro = items[4] + if nro != 0: + params, scattering_radius = get_tab1_record(file_obj) + + # Get SPI, AP, and LSSF + formalism = items[3] + if not (fission_widths and formalism == 1): + items = get_cont_record(file_obj) + target_spin = items[0] + if nro == 0: + scattering_radius = items[1] + add_to_background = (items[2] == 0) + + if not fission_widths and formalism == 1: + # Case A -- fission widths not given, all parameters are + # energy-independent + NLS = items[4] + columns = ['L', 'J', 'd', 'amun', 'gn', 'gg', 'gf'] + records = [] + for ls in range(NLS): + items, values = get_list_record(file_obj) + l = items[2] + NJS = items[5] + for j in range(NJS): + d, j, amun, gn, gg = values[6*j:6*j + 5] + records.append([l, j, d, amun, gn, gg, 0.0]) + parameters = pd.DataFrame.from_records(records, columns=columns) + energies = None + + elif fission_widths and formalism == 1: + # Case B -- fission widths given, only fission widths are + # energy-dependent + items, energies = get_list_record(file_obj) + target_spin = items[0] + if nro == 0: + scattering_radius = items[1] + add_to_background = (items[2] == 0) + NE, NLS = items[4:6] + l_values = np.zeros(NLS, int) + records = [[] for e in energies] + columns = ['L', 'J', 'd', 'amun', 'gn0', 'gg', 'gf'] + for ls in range(NLS): + items = get_cont_record(file_obj) + l = items[2] + NJS = items[4] + for j in range(NJS): + items, values = get_list_record(file_obj) + muf = items[3] + d = values[0] + j = values[1] + amun = values[2] + gn0 = values[3] + gg = values[4] + gf = values[6:] + for k, gf_i in enumerate(gf): + records[k].append([l, j, d, amun, gn0, gg, gf]) + parameters = [pd.DataFrame(r, columns=columns) for r in records] + + elif formalism == 2: + # Case C -- all parameters are energy-dependent + NLS = items[4] + columns = ['L', 'J', 'E', 'd', 'gx', 'gn0', 'gg', 'gf'] + records = [] + for ls in range(NLS): + items = get_cont_record(file_obj) + l = items[2] + NJS = items[4] + for j in range(NJS): + items, values = get_list_record(file_obj) + ne = items[5] + j = items[0] + amux = values[2] + amun = values[3] + amug = values[4] + amuf = values[5] + for k in range(1, ne + 1): + E = values[6*k] + d = values[6*k + 1] + gx = values[6*k + 2] + gn0 = values[6*k + 3] + gg = values[6*k + 4] + gf = values[6*k + 5] + records.append([l, j, E, d, gx, gn0, gg, gf]) + parameters = pd.DataFrame.from_records(records, columns=columns) + energies = None + + urr = cls(target_spin, energy_min, energy_max, scattering_radius) + urr.parameters = parameters + urr.add_to_background = add_to_background + urr.energies = energies + + return urr + + +_FORMALISMS = {0: ResonanceRange, + 1: SingleLevelBreitWigner, + 2: MultiLevelBreitWigner, + 3: ReichMoore, + 7: RMatrixLimited} + +_RESOLVED = (SingleLevelBreitWigner, MultiLevelBreitWigner, + ReichMoore, RMatrixLimited) diff --git a/setup.py b/setup.py index 4479bb82c..52ffaecf4 100644 --- a/setup.py +++ b/setup.py @@ -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)