From cbcd3a7f573c1c5ead8a7d1ec5457d34fa034f87 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 24 Feb 2016 11:01:12 -0600 Subject: [PATCH 01/32] 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) From eed74c9d2632944f041196a19f0177088de5637a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Sep 2016 14:17:25 -0500 Subject: [PATCH 02/32] Fix issue with losing precision in RM reconstruction when K and phi are small --- openmc/data/reconstruct.pyx | 47 +++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx index 4790cc345..452968d14 100644 --- a/openmc/data/reconstruct.pyx +++ b/openmc/data/reconstruct.pyx @@ -66,7 +66,6 @@ cdef double phaseshift(int l, double rho): Hardsphere phase shift """ - if l == 0: return rho elif l == 1: @@ -377,7 +376,8 @@ def reconstruct_rm(rm, double E): 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 double Kr, Ki, x + cdef double complex Ubar, U_, factor cdef bint hasfission cdef np.ndarray[double, ndim=2] one cdef np.ndarray[double complex, ndim=2] K, Imat, U @@ -459,31 +459,24 @@ def reconstruct_rm(rm, double E): # 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 + # Calculate j/2 * inverse of denominator of K matrix terms + factor = 0.5j/(E_r - E - 0.5j*gg) # Upper triangular portion of K matrix -- see ENDF-102, # Equation D.28 - K[0,0] = K[0,0] + gn*gn*denominator_inverse + K[0,0] = K[0,0] + gn*gn*factor 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 + K[0,1] = K[0,1] + gn*gfa*factor + K[0,2] = K[0,2] + gn*gfb*factor + K[1,1] = K[1,1] + gfa*gfa*factor + K[1,2] = K[1,2] + gfa*gfb*factor + K[2,2] = K[2,2] + gfb*gfb*factor 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: @@ -500,9 +493,23 @@ def reconstruct_rm(rm, double E): # 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 + U_ = Ubar*(2/(1 - K[0,0]) - 1) + if abs(creal(K[0,0])) < 3e-4 and abs(phi) < 3e-4: + # If K and phi are both very small, the calculated cross + # sections can lose precision because the real part of U + # ends up very close to unity. To get around this, we + # use Euler's formula to express Ubar by real and + # imaginary parts, expand cos(2phi) = 1 - 2phi^2 + + # O(phi^4), and then simplify + Kr = creal(K[0,0]) + Ki = cimag(K[0,0]) + x = 2*(-Kr + (Kr*Kr + Ki*Ki)*(1 - phi*phi) + phi*phi - + sin(2*phi)*Ki)/((1 - Kr)*(1 - Kr) + Ki*Ki) + total += 2*gJ*x + elastic += gJ*(x*x + cimag(U_)**2) + else: + total += 2*gJ*(1 - creal(U_)) # ENDF-102, Eq. D.23 + elastic += gJ*cabs(1 - U_)**2 # ENDF-102, Eq. D.24 # Calculate capture as difference of other cross sections as per ENDF-102, # Equation D.25 From f08550f5ed830fbe39786b490c94a2b3b8618705 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Sep 2016 14:32:08 -0500 Subject: [PATCH 03/32] Handle reconstruction for multiple resolved ranges (Pu239 in ENDF/B-VII.1) --- openmc/data/function.py | 25 ++++++++++++++----------- openmc/data/neutron.py | 2 +- openmc/data/resonance.py | 4 ++++ 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 13df12b47..4121918e2 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -459,18 +459,21 @@ class ResonancesWithBackground(EqualityMixin): # 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) + for r in self.resonances: + if not isinstance(r, openmc.data.resonance._RESOLVED): + continue - # 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] + if isinstance(x, Iterable): + # Determine which energies are within resolved resonance range + within = (r.energy_min <= x) & (x <= r.energy_max) + + # Get resonance cross sections and add to background + resonant_xs = r.reconstruct(x[within]) + xs[within] += resonant_xs[self.mt] + else: + if r.energy_min <= x <= r.energy_max: + resonant_xs = r.reconstruct(x) + xs += resonant_xs[self.mt] return xs diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 0c74aeaa1..f903ea6f5 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -659,7 +659,7 @@ class IncidentNeutron(EqualityMixin): # Replace cross sections for elastic, capture, fission try: - if isinstance(data.resonances.resolved, _RESOLVED): + if any(isinstance(r, _RESOLVED) for r in data.resonances): for mt in (2, 102, 18): if mt in data.reactions: rx = data.reactions[mt] diff --git a/openmc/data/resonance.py b/openmc/data/resonance.py index 206a56df4..77d522a04 100644 --- a/openmc/data/resonance.py +++ b/openmc/data/resonance.py @@ -36,6 +36,10 @@ class Resonances(object): def __init__(self, ranges): self.ranges = ranges + def __iter__(self): + for r in self.ranges: + yield r + @property def ranges(self): return self._ranges From 5555cfd1aaefcccdc9dcaeba847dc8d38ee21dc0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 30 Sep 2016 14:44:44 -0500 Subject: [PATCH 04/32] Respond to @smharper comments on #721 --- openmc/data/nbody.py | 4 ++-- openmc/data/reaction.py | 4 +--- openmc/data/reconstruct.pyx | 4 ++-- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/openmc/data/nbody.py b/openmc/data/nbody.py index 6616b1858..b186c253a 100644 --- a/openmc/data/nbody.py +++ b/openmc/data/nbody.py @@ -18,7 +18,7 @@ class NBodyPhaseSpace(AngleEnergy): atomic_weight_ratio : float Atomic weight ratio of target nuclide q_value : float - Q value for reaction in MeV + Q value for reaction in MeV or eV, depending on the data source. Attributes ---------- @@ -29,7 +29,7 @@ class NBodyPhaseSpace(AngleEnergy): atomic_weight_ratio : float Atomic weight ratio of target nuclide q_value : float - Q value for reaction in MeV + Q value for reaction in MeV or eV, depending on the data source. """ diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index cef635135..fa55fbab9 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -663,9 +663,7 @@ class Reaction(EqualityMixin): mt : int The ENDF MT number for this reaction. q_value : float - The Q-value of this reaction in MeV. - threshold : float - Threshold of the reaction in MeV + The Q-value of this reaction in MeV or eV, depending on the data source. xs : dict of str to openmc.data.Function1D Microscopic cross section for this reaction as a function of incident energy; these cross sections are provided in a dictionary where the key diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx index 452968d14..f63a155b1 100644 --- a/openmc/data/reconstruct.pyx +++ b/openmc/data/reconstruct.pyx @@ -15,8 +15,8 @@ cdef extern from "complex.h": 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 +cdef double NEUTRON_MASS_ENERGY = 939.5654133e6 # eV/c^2 +cdef double HBAR_C = 197.3269788e5 # eV-b^0.5 @cython.cdivision(True) From 49fb1ade00eb0c6c5c0b842f66c7622f37b715bc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 2 Oct 2016 15:59:57 -0400 Subject: [PATCH 05/32] Allow wmp to work with s(a, b) and URR interp --- docs/source/usersguide/input.rst | 16 ++++++++-------- openmc/settings.py | 24 ++++++++++++++++++++++-- src/constants.F90 | 3 +-- src/cross_section.F90 | 2 -- src/global.F90 | 3 +++ src/input_xml.F90 | 21 ++++++++++++++++----- src/nuclide_header.F90 | 4 ---- src/physics.F90 | 2 +- src/sab_header.F90 | 4 ---- tests/test_multipole/inputs_true.dat | 2 +- tests/test_multipole/results_true.dat | 2 +- tests/test_multipole/test_multipole.py | 4 +++- 12 files changed, 56 insertions(+), 31 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 9b2d78607..2e4d3f087 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -292,8 +292,8 @@ OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range cross sections. If this element is absent from the settings.xml file, the :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. - .. note:: The :ref:`temperature_method` must also be set to "multipole" for - windowed multipole functionality. + .. note:: The element must also be set to "true" + for windowed multipole functionality. ```` Element --------------------------- @@ -725,16 +725,16 @@ a material default temperature. ```` Element -------------------------------- -The ```` element has an accepted value of "nearest", -"interpolation", or "multipole". A value of "nearest" indicates that for each +The ```` element has an accepted value of "nearest" or +"interpolation". A value of "nearest" indicates that for each cell, the nearest temperature at which cross sections are given is to be applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of "interpolation" indicates that cross sections are to be linear-linear interpolated between temperatures at which nuclear data are present (see -:ref:`temperature_treatment`). A value of "multipole" indicates that the -windowed multipole method should be used to evaluate temperature-dependent cross -sections in the resolved resonance range (a :ref:`windowed multipole library -` must also be available). +:ref:`temperature_treatment`). Note that if ```` is set +true, then for which multipole data is available will use the windowed multipole +method for cross sections in the resolved resonance range. (a +:ref:`windowed multipole library ` must also be available). *Default*: "nearest" diff --git a/openmc/settings.py b/openmc/settings.py index f255a8bb1..0d3eb61db 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -105,7 +105,7 @@ class Settings(object): temperatures at which nuclear data doesn't exist. Accepted keys are 'default', 'method', and 'tolerance'. The value for 'default' should be a float representing the default temperature in Kelvin. The value for - 'method' should be 'nearest' or 'multipole'. If the method is + 'method' should be 'nearest' or 'interpolation'. If the method is 'nearest', 'tolerance' indicates a range of temperature within which cross sections may be used. trigger_active : bool @@ -135,6 +135,9 @@ class Settings(object): Coordinates of the lower-left point of the UFS mesh ufs_upper_right : tuple or list Coordinates of the upper-right point of the UFS mesh + use_windowed_multipole : bool + Whether or not windowed multipole can be used to evaluate resolved + resonance cross sections. resonance_scattering : ResonanceScattering or iterable of ResonanceScattering The elastic scattering model to use for resonant isotopes volume_calculations : VolumeCalculation or iterable of VolumeCalculation @@ -219,6 +222,7 @@ class Settings(object): self._settings_file = ET.Element("settings") self._run_mode_subelement = None + self._multipole_active = None self._resonance_scattering = cv.CheckedList( ResonanceScattering, 'resonance scattering models') @@ -417,6 +421,10 @@ class Settings(object): def dd_count_interactions(self): return self._dd_count_interactions + @property + def use_windowed_multipole(self): + return self._multipole_active + @property def resonance_scattering(self): return self._resonance_scattering @@ -676,7 +684,7 @@ class Settings(object): cv.check_type('default temperature', value, Real) elif key == 'method': cv.check_value('temperature method', value, - ['nearest', 'interpolation', 'multipole']) + ['nearest', 'interpolation']) elif key == 'tolerance': cv.check_type('temperature tolerance', value, Real) self._temperature = temperature @@ -808,6 +816,11 @@ class Settings(object): self._dd_count_interactions = interactions + @use_windowed_multipole.setter + def use_windowed_multipole(self, active): + cv.check_type('use_windowed_multipole', active, bool) + self._multipole_active = active + @resonance_scattering.setter def resonance_scattering(self, res): if not isinstance(res, MutableSequence): @@ -1111,6 +1124,12 @@ class Settings(object): subelement = ET.SubElement(element, "count_interactions") subelement.text = str(self._dd_count_interactions).lower() + def _create_use_multipole_subelement(self): + if self._multipole_active is not None: + element = ET.SubElement(self._settings_file, + "use_windowed_multipole") + element.text = str(self._multipole_active) + def _create_resonance_scattering_subelement(self): if len(self.resonance_scattering) > 0: elem = ET.SubElement(self._settings_file, 'resonance_scattering') @@ -1158,6 +1177,7 @@ class Settings(object): self._create_track_subelement() self._create_ufs_subelement() self._create_dd_subelement() + self._create_use_multipole_subelement() self._create_resonance_scattering_subelement() self._create_volume_calcs_subelement() diff --git a/src/constants.F90 b/src/constants.F90 index cc0d663ea..8c321156c 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -271,8 +271,7 @@ module constants ! Temperature treatment method integer, parameter :: & TEMPERATURE_NEAREST = 1, & - TEMPERATURE_INTERPOLATION = 2, & - TEMPERATURE_MULTIPOLE = 3 + TEMPERATURE_INTERPOLATION = 2 ! ============================================================================ ! TALLY-RELATED CONSTANTS diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 4bb323c51..69363e40a 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -185,8 +185,6 @@ contains f = (kT - nuc % kTs(i_temp)) / & (nuc % kTs(i_temp + 1) - nuc % kTs(i_temp)) if (f > prn()) i_temp = i_temp + 1 - case (TEMPERATURE_MULTIPOLE) - i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) end select end if diff --git a/src/global.F90 b/src/global.F90 index bea5f61a8..d424e1cec 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -99,6 +99,9 @@ module global ! What to assume for expanding natural elements integer :: default_expand = ENDF_BVII1 + ! Whether or not windowed multipole cross sections should be used. + logical :: multipole_active = .false. + ! Default temperature and method for choosing temperatures integer :: temperature_method = TEMPERATURE_NEAREST real(8) :: temperature_tolerance = 10.0_8 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 16892397a..74b5f47ae 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1063,6 +1063,20 @@ contains end select end if + ! Check to see if windowed multipole functionality is requested + if (check_for_node(doc, "use_windowed_multipole")) then + call get_node_value(doc, "use_windowed_multipole", temp_str) + select case (to_lower(temp_str)) + case ('true', '1') + multipole_active = .true. + case ('false', '0') + multipole_active = .false. + case default + call fatal_error("Unrecognized value for in & + &settings.xml") + end select + end if + call get_node_list(doc, "volume_calc", node_vol_list) n = get_list_size(node_vol_list) allocate(volume_calcs(n)) @@ -1082,8 +1096,6 @@ contains temperature_method = TEMPERATURE_NEAREST case ('interpolation') temperature_method = TEMPERATURE_INTERPOLATION - case ('multipole') - temperature_method = TEMPERATURE_MULTIPOLE case default call fatal_error("Unknown temperature method: " // trim(temp_str)) end select @@ -5786,8 +5798,7 @@ contains call already_read % add(name) ! Read multipole file into the appropriate entry on the nuclides array - if (temperature_method == TEMPERATURE_MULTIPOLE) & - call read_multipole_data(i_nuclide) + if (multipole_active) call read_multipole_data(i_nuclide) end if ! Check if material is fissionable @@ -5841,7 +5852,7 @@ contains end do ! If the user wants multipole, make sure we found a multipole library. - if (temperature_method == TEMPERATURE_MULTIPOLE) then + if (multipole_active) then mp_found = .false. do i = 1, size(nuclides) if (nuclides(i) % mp_present) then diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 4235f4234..d07ccbe5a 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -302,10 +302,6 @@ module nuclide_header trim(to_str(nint(temp_desired))) // " K.") end do TEMP_LOOP - case (TEMPERATURE_MULTIPOLE) - ! Add first available temperature - call temps_to_read % push_back(nint(temps_available(1))) - end select ! Sort temperatures to read diff --git a/src/physics.F90 b/src/physics.F90 index 23fb6f1df..6a6a9787d 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -334,7 +334,7 @@ contains else ! Determine temperature - if (temperature_method == TEMPERATURE_MULTIPOLE) then + if (nuc % mp_present) then kT = p % sqrtkT**2 else kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 21a2ab032..dcf5df979 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -187,10 +187,6 @@ contains trim(to_str(nint(temp_desired))) // " K.") end do TEMP_LOOP - case (TEMPERATURE_MULTIPOLE) - ! Add first available temperature - call temps_to_read % push_back(nint(temps_available(1))) - end select ! Sort temperatures to read diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 930be9536..9ff8c6e39 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -1 +1 @@ -8462e17d102259b3a48a7e908bc75038a28a19d4a5e8bd38f26579e5baf3998bc2d05d1d3055ac1ed478d0c01bd64868d654919c2e6ca3fea86d79f03eda25ef \ No newline at end of file +cab3356c163ceace74251d20cafc1b46f9fc3af428685cf5cea9964f5356d59a967468a14b3a2cd2ad3eea04b1d164b5daee8259f38be8ec5fb04fa7b4982830 \ No newline at end of file diff --git a/tests/test_multipole/results_true.dat b/tests/test_multipole/results_true.dat index 4d379b3f4..b17f6b0aa 100644 --- a/tests/test_multipole/results_true.dat +++ b/tests/test_multipole/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.425673E+00 1.779969E-02 +1.363786E+00 1.103929E-02 Cell ID = 11 Name = diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index a44489aa5..e908eba2d 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -17,6 +17,7 @@ class MultipoleTestHarness(PyAPITestHarness): moderator.set_density('g/cc', 1.0) moderator.add_nuclide('H1', 2.0) moderator.add_nuclide('O16', 1.0) + moderator.add_s_alpha_beta('c_H_in_H2O') dense_fuel = openmc.Material(material_id=2) dense_fuel.set_density('g/cc', 4.5) @@ -67,7 +68,8 @@ class MultipoleTestHarness(PyAPITestHarness): sets_file.particles = 1000 sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1])) sets_file.output = {'summary': True} - sets_file.temperature = {'method': 'multipole'} + sets_file.use_windowed_multipole = True + sets_file.temperature = {'tolerance': 1000} sets_file.export_to_xml() #################### From a7d8a49e07c2ed3aa9eba94044d4b92762be9611 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Oct 2016 23:08:34 -0500 Subject: [PATCH 06/32] Describe skipC argument in get_cont_record --- openmc/data/endf.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index ba398598d..0940fcc8e 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -104,6 +104,9 @@ def get_cont_record(file_obj, skipC=False): ---------- file_obj : file-like object ENDF-6 file to read from + skipC : bool + Determine whether to skip the first two quantities (C1, C2) of the CONT + record. Returns ------- From cbd2637efc4dfdd52ce3f4dc6469dba761cdcbaf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Oct 2016 15:08:59 -0500 Subject: [PATCH 07/32] Handle energy-dependent delayed neutron group probabilities in ENDF files. --- openmc/data/reaction.py | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index fa55fbab9..081cbc245 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -419,8 +419,37 @@ def _get_fission_products_endf(ev): dist = UncorrelatedAngleEnergy() dist.energy = EnergyDistribution.from_endf(file_obj, params) - delayed_neutron = products[-6 + i] - delayed_neutron.yield_.y *= applicability.y[0] + delayed_neutron = products[1 + i] + yield_ = delayed_neutron.yield_ + + # Here we handle the fact that the delayed neutron yield is the + # product of the total delayed neutron yield and the + # "applicability" of the energy distribution law in file 5. + if isinstance(yield_, Tabulated1D): + if np.all(applicability.y == applicability.y[0]): + yield_.y *= applicability.y[0] + else: + # Get union energy grid and ensure energies are within + # interpolable range of both functions + max_energy = min(yield_.x[-1], applicability.x[-1]) + energy = np.union1d(yield_.x, applicability.x) + energy = energy[energy <= max_energy] + + # Calculate group yield + group_yield = yield_(energy) * applicability(energy) + delayed_neutron.yield_ = Tabulated1D(energy, group_yield) + elif isinstance(yield_, Polynomial): + if len(yield_) == 1: + delayed_neutron.yield_ = deepcopy(applicability) + delayed_neutron.yield_.y *= yield_.coef[0] + else: + if np.all(applicability.y == applicability.y[0]): + yield_.coef[0] *= applicability.y[0] + else: + raise NotImplementedError( + 'Total delayed neutron yield and delayed group ' + 'probability are both energy-dependent.') + delayed_neutron.distribution.append(dist) return products, derived_products From b3e351474c50f7cb5283890a5f2c5fa912c06391 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Oct 2016 15:10:12 -0500 Subject: [PATCH 08/32] Handle Sher-Beck energy release data with polynomial nu. --- openmc/data/fission_energy.py | 38 +++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 15 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index c7f71b657..9479e6731 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -405,27 +405,35 @@ class FissionEnergyRelease(EqualityMixin): # MT=18 (n, fission) might not be available so try MT=19 (n, f) as # well. if 18 in incident_neutron.reactions: - nu_prompt = [p for p in incident_neutron[18].products - if p.particle == 'neutron' - and p.emission_mode == 'prompt'] + nu = [p.yield_ for p in incident_neutron[18].products + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] elif 19 in incident_neutron.reactions: - nu_prompt = [p for p in incident_neutron[19].products - if p.particle == 'neutron' - and p.emission_mode == 'prompt'] + nu = [p.yield_ for p in incident_neutron[19].products + if p.particle == 'neutron' + and p.emission_mode in ('prompt', 'total')] else: raise ValueError('IncidentNeutron data has no fission ' 'reaction.') - if len(nu_prompt) == 0: + if len(nu) == 0: raise ValueError('Nu data is needed to compute fission energy ' 'release with the Sher-Beck format.') - if len(nu_prompt) > 1: - raise ValueError('Ambiguous prompt value.') - if not isinstance(nu_prompt[0].yield_, Tabulated1D): - raise TypeError('Sher-Beck fission energy release currently ' - 'only supports Tabulated1D nu data.') - ENP = deepcopy(nu_prompt[0].yield_) - ENP.y = (energy_release['ENP'] + 1.307 * ENP.x - - 8.07 * (ENP.y - ENP.y[0])) + if len(nu) > 1: + raise ValueError('Ambiguous prompt/total nu value.') + + nu = nu[0] + if isinstance(nu, Tabulated1D): + ENP = deepcopy(nu) + ENP.y = (energy_release['ENP'] + 1.307 * nu.x + - 8.07 * (nu.y - nu.y[0])) + elif isinstance(nu, Polynomial): + if len(nu) == 1: + ENP = Polynomial([energy_release['ENP'][0], 1.307]) + else: + ENP = Polynomial( + [energy_release['ENP'][0], 1.307 - 8.07*nu.coef[1]] + + [-8.07*c for c in nu.coef[2:]]) + out.prompt_neutrons = ENP return out From 9079ae7b7478919d2cba685d3cdf491f6bc13738 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Oct 2016 15:10:41 -0500 Subject: [PATCH 09/32] Ignore obsolete energy transformation matrix in MF=4. --- openmc/data/angle_distribution.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/openmc/data/angle_distribution.py b/openmc/data/angle_distribution.py index 0f0105c2a..669951059 100644 --- a/openmc/data/angle_distribution.py +++ b/openmc/data/angle_distribution.py @@ -1,6 +1,7 @@ from collections import Iterable from io import StringIO from numbers import Real +from warnings import warn import numpy as np @@ -224,13 +225,23 @@ class AngleDistribution(EqualityMixin): # Read HEAD record items = get_head_record(file_obj) + lvt = items[2] ltt = items[3] # Read CONT record items = get_cont_record(file_obj) li = items[2] + nk = items[4] center_of_mass = (items[3] == 2) + # Check for obsolete energy transformation matrix. If present, just skip + # it and keep reading + if lvt > 0: + warn('Obsolete energy transformation matrix in MF=4 angular ' + 'distribution.') + for _ in range((nk + 5)//6): + file_obj.readline() + if ltt == 0 and li == 1: # Purely isotropic energy = np.array([0., ev.info['energy_max']]) From 348685057a3638bed5b64faa392f43da5ad9b717 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Oct 2016 16:36:55 -0500 Subject: [PATCH 10/32] Handle missing MOD values in JEFF 3.2 gracefully --- openmc/data/endf.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index 0940fcc8e..0c590be0a 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -384,10 +384,17 @@ class Evaluation(object): # 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)) - + line = file_obj.readline() + mf = int(line[22:33]) + mt = int(line[33:44]) + nc = int(line[44:55]) + try: + mod = int(line[55:66]) + except ValueError: + # In JEFF 3.2, a few isotopes of U have MOD values that are + # missing. This prevents failure on these isotopes. + mod = 0 + self.reaction_list.append((mf, mt, nc, mod)) class Tabulated2D(object): From e2ff3a83572a00383f8e56aa629a0ecee13b1c22 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Oct 2016 16:47:33 -0500 Subject: [PATCH 11/32] Add error message when number of delayed fission spectra doesn't match number of precursor groups. --- openmc/data/reaction.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 081cbc245..af903fe8a 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -413,7 +413,11 @@ def _get_fission_products_endf(ev): file_obj = StringIO(ev.section[5, 455]) items = get_head_record(file_obj) nk = items[4] - assert nk == len(decay_constants) + if nk != len(decay_constants): + raise ValueError( + 'Number of delayed neutron fission spectra ({}) does not ' + 'match number of delayed neutron precursors ({}).'.format( + nk, len(decay_constants))) for i in range(nk): params, applicability = get_tab1_record(file_obj) dist = UncorrelatedAngleEnergy() From 36cf8934fd5ed9875963439eddd0bea76557dc74 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Oct 2016 22:17:52 -0500 Subject: [PATCH 12/32] Fix units for fission energy release in ENDF files --- openmc/data/fission_energy.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 9479e6731..687cb5703 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -16,13 +16,15 @@ if sys.version_info[0] >= 3: basestring = str -def _extract_458_data(ev): +def _extract_458_data(ev, units='eV'): """Read an ENDF file and extract the MF=1, MT=458 values. Parameters ---------- ev : openmc.data.Evaluation ENDF evaluation + units : {'eV', 'MeV'} + The units are used in values returned. Returns ------- @@ -86,12 +88,13 @@ def _extract_458_data(ev): for coeffs in uncertainty.values(): coeffs[2] *= 1e-6 # Convert eV to MeV. - for coeffs in value.values(): - for i in range(len(coeffs)): - coeffs[i] *= 10**(-6 + 6*i) - for coeffs in uncertainty.values(): - for i in range(len(coeffs)): - coeffs[i] *= 10**(-6 + 6*i) + if units == 'MeV': + for coeffs in value.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) + for coeffs in uncertainty.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) return value, uncertainty @@ -158,7 +161,7 @@ def write_compact_458_library(endf_files, output_name='fission_Q_data.h5', continue # Get the important bits. - data = _extract_458_data(ev) + data = _extract_458_data(ev, 'MeV') if data is None: continue value, uncertainty = data From b68e3c12aca9af941b624116e22d2522342af05c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Oct 2016 13:52:07 -0400 Subject: [PATCH 13/32] Stick multipole flag back in settings.temperature --- openmc/settings.py | 44 ++++++++++---------------- tests/test_multipole/test_multipole.py | 3 +- 2 files changed, 17 insertions(+), 30 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 0d3eb61db..5019980eb 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -103,11 +103,13 @@ class Settings(object): temperature : dict Defines a default temperature and method for treating intermediate temperatures at which nuclear data doesn't exist. Accepted keys are - 'default', 'method', and 'tolerance'. The value for 'default' should be - a float representing the default temperature in Kelvin. The value for - 'method' should be 'nearest' or 'interpolation'. If the method is - 'nearest', 'tolerance' indicates a range of temperature within which - cross sections may be used. + 'default', 'method', 'tolerance', and 'multipole'. The value for + 'default' should be a float representing the default temperature in + Kelvin. The value for 'method' should be 'nearest' or 'interpolation'. + If the method is 'nearest', 'tolerance' indicates a range of temperature + within which cross sections may be used. 'multipole' is a boolean + indicating whether or not the windowed multipole method should be used + to evaluate resolved resonance cross sections. trigger_active : bool Indicate whether tally triggers are used trigger_max_batches : int @@ -135,9 +137,6 @@ class Settings(object): Coordinates of the lower-left point of the UFS mesh ufs_upper_right : tuple or list Coordinates of the upper-right point of the UFS mesh - use_windowed_multipole : bool - Whether or not windowed multipole can be used to evaluate resolved - resonance cross sections. resonance_scattering : ResonanceScattering or iterable of ResonanceScattering The elastic scattering model to use for resonant isotopes volume_calculations : VolumeCalculation or iterable of VolumeCalculation @@ -222,7 +221,6 @@ class Settings(object): self._settings_file = ET.Element("settings") self._run_mode_subelement = None - self._multipole_active = None self._resonance_scattering = cv.CheckedList( ResonanceScattering, 'resonance scattering models') @@ -421,10 +419,6 @@ class Settings(object): def dd_count_interactions(self): return self._dd_count_interactions - @property - def use_windowed_multipole(self): - return self._multipole_active - @property def resonance_scattering(self): return self._resonance_scattering @@ -679,7 +673,7 @@ class Settings(object): cv.check_type('temperature settings', temperature, Mapping) for key, value in temperature.items(): cv.check_value('temperature key', key, - ['default', 'method', 'tolerance']) + ['default', 'method', 'tolerance', 'multipole']) if key == 'default': cv.check_type('default temperature', value, Real) elif key == 'method': @@ -687,6 +681,8 @@ class Settings(object): ['nearest', 'interpolation']) elif key == 'tolerance': cv.check_type('temperature tolerance', value, Real) + elif key == 'multipole': + cv.check_type('temperature multipole', value, bool) self._temperature = temperature @threads.setter @@ -816,11 +812,6 @@ class Settings(object): self._dd_count_interactions = interactions - @use_windowed_multipole.setter - def use_windowed_multipole(self, active): - cv.check_type('use_windowed_multipole', active, bool) - self._multipole_active = active - @resonance_scattering.setter def resonance_scattering(self, res): if not isinstance(res, MutableSequence): @@ -1063,8 +1054,12 @@ class Settings(object): def _create_temperature_subelements(self): if self.temperature: for key, value in self.temperature.items(): - element = ET.SubElement(self._settings_file, - "temperature_{}".format(key)) + if (key == 'multipole'): + element = ET.SubElement(self._settings_file, + "use_windowed_multipole") + else: + element = ET.SubElement(self._settings_file, + "temperature_{}".format(key)) element.text = str(value) def _create_threads_subelement(self): @@ -1124,12 +1119,6 @@ class Settings(object): subelement = ET.SubElement(element, "count_interactions") subelement.text = str(self._dd_count_interactions).lower() - def _create_use_multipole_subelement(self): - if self._multipole_active is not None: - element = ET.SubElement(self._settings_file, - "use_windowed_multipole") - element.text = str(self._multipole_active) - def _create_resonance_scattering_subelement(self): if len(self.resonance_scattering) > 0: elem = ET.SubElement(self._settings_file, 'resonance_scattering') @@ -1177,7 +1166,6 @@ class Settings(object): self._create_track_subelement() self._create_ufs_subelement() self._create_dd_subelement() - self._create_use_multipole_subelement() self._create_resonance_scattering_subelement() self._create_volume_calcs_subelement() diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index e908eba2d..32b38dd76 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -68,8 +68,7 @@ class MultipoleTestHarness(PyAPITestHarness): sets_file.particles = 1000 sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1])) sets_file.output = {'summary': True} - sets_file.use_windowed_multipole = True - sets_file.temperature = {'tolerance': 1000} + sets_file.temperature = {'tolerance': 1000, 'multipole': True} sets_file.export_to_xml() #################### From 63b9543fc9aa93d2c0302c2b8c9827f3c6ce0706 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Oct 2016 14:19:17 -0400 Subject: [PATCH 14/32] use_windowed_multipole --> temperature_multipole --- docs/source/usersguide/input.rst | 33 ++++++++++++++-------------- openmc/settings.py | 8 ++----- src/global.F90 | 4 +--- src/input_xml.F90 | 30 ++++++++++++------------- src/relaxng/settings.rnc | 4 ++-- src/relaxng/settings.rng | 10 ++++----- tests/test_multipole/inputs_true.dat | 2 +- 7 files changed, 41 insertions(+), 50 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 2e4d3f087..847dc66ff 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -292,8 +292,8 @@ OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range cross sections. If this element is absent from the settings.xml file, the :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. - .. note:: The element must also be set to "true" - for windowed multipole functionality. + .. note:: The element must also be set to "true" for + windowed multipole functionality. ```` Element --------------------------- @@ -731,13 +731,23 @@ cell, the nearest temperature at which cross sections are given is to be applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of "interpolation" indicates that cross sections are to be linear-linear interpolated between temperatures at which nuclear data are present (see -:ref:`temperature_treatment`). Note that if ```` is set -true, then for which multipole data is available will use the windowed multipole -method for cross sections in the resolved resonance range. (a -:ref:`windowed multipole library ` must also be available). +:ref:`temperature_treatment`). *Default*: "nearest" +.. _temperature_multipole: + +```` Element +----------------------------------- + +The ```` element toggles the windowed multipole +capability on or off. If this element is set to "True" and the relevant data is +available, OpenMC will use the windowed multipole method to evaluate and Doppler +broaden cross sections in the resolved resonance range. This override other +methods like "nearest" and "interpolation" in the resolved resonance range. + + *Default*: False + .. _temperature_tolerance: ```` Element @@ -850,17 +860,6 @@ problem. It has the following attributes/sub-elements: *Default*: None - -```` Element ------------------------------------- - -The ```` element toggles the windowed multipole -capability on or off. If this element is set to "True" and the relevant data is -available, OpenMC will use the windowed multipole method to evaluate and Doppler -broaden cross sections in the resolved resonance range. - - *Default*: False - ```` Element ----------------------- diff --git a/openmc/settings.py b/openmc/settings.py index 5019980eb..84fe20485 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1054,12 +1054,8 @@ class Settings(object): def _create_temperature_subelements(self): if self.temperature: for key, value in self.temperature.items(): - if (key == 'multipole'): - element = ET.SubElement(self._settings_file, - "use_windowed_multipole") - else: - element = ET.SubElement(self._settings_file, - "temperature_{}".format(key)) + element = ET.SubElement(self._settings_file, + "temperature_{}".format(key)) element.text = str(value) def _create_threads_subelement(self): diff --git a/src/global.F90 b/src/global.F90 index d424e1cec..c41ab52c0 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -99,11 +99,9 @@ module global ! What to assume for expanding natural elements integer :: default_expand = ENDF_BVII1 - ! Whether or not windowed multipole cross sections should be used. - logical :: multipole_active = .false. - ! Default temperature and method for choosing temperatures integer :: temperature_method = TEMPERATURE_NEAREST + logical :: temperature_multipole = .false. real(8) :: temperature_tolerance = 10.0_8 real(8) :: temperature_default = 293.6_8 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 74b5f47ae..8691a9d5c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1063,20 +1063,6 @@ contains end select end if - ! Check to see if windowed multipole functionality is requested - if (check_for_node(doc, "use_windowed_multipole")) then - call get_node_value(doc, "use_windowed_multipole", temp_str) - select case (to_lower(temp_str)) - case ('true', '1') - multipole_active = .true. - case ('false', '0') - multipole_active = .false. - case default - call fatal_error("Unrecognized value for in & - &settings.xml") - end select - end if - call get_node_list(doc, "volume_calc", node_vol_list) n = get_list_size(node_vol_list) allocate(volume_calcs(n)) @@ -1103,6 +1089,18 @@ contains if (check_for_node(doc, "temperature_tolerance")) then call get_node_value(doc, "temperature_tolerance", temperature_tolerance) end if + if (check_for_node(doc, "temperature_multipole")) then + call get_node_value(doc, "temperature_multipole", temp_str) + select case (to_lower(temp_str)) + case ('true', '1') + temperature_multipole = .true. + case ('false', '0') + temperature_multipole = .false. + case default + call fatal_error("Unrecognized value for in & + &settings.xml") + end select + end if ! Close settings XML file call close_xmldoc(doc) @@ -5798,7 +5796,7 @@ contains call already_read % add(name) ! Read multipole file into the appropriate entry on the nuclides array - if (multipole_active) call read_multipole_data(i_nuclide) + if (temperature_multipole) call read_multipole_data(i_nuclide) end if ! Check if material is fissionable @@ -5852,7 +5850,7 @@ contains end do ! If the user wants multipole, make sure we found a multipole library. - if (multipole_active) then + if (temperature_multipole) then mp_found = .false. do i = 1, size(nuclides) if (nuclides(i) % mp_present) then diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 70550a40f..b05ab7d6b 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -132,6 +132,8 @@ element settings { element temperature_method { xsd:string }? & + element temperature_multipole { xsd:boolean }? & + element temperature_tolerance { xsd:double }? & element threads { xsd:positiveInteger }? & @@ -182,6 +184,4 @@ element settings { attribute E_max { xsd:double })? }* }? & - - element use_windowed_multipole { xsd:boolean }? } diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 246c78e68..cbeac25a5 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -575,6 +575,11 @@ + + + + + @@ -816,10 +821,5 @@ - - - - - diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 9ff8c6e39..3b2b43ec7 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -1 +1 @@ -cab3356c163ceace74251d20cafc1b46f9fc3af428685cf5cea9964f5356d59a967468a14b3a2cd2ad3eea04b1d164b5daee8259f38be8ec5fb04fa7b4982830 \ No newline at end of file +2be927608035759f52a2ac88b2c84e28c06e0f7905187bfd4695fecd80f417e727d8879c52fe03253938f7aa0301061f72dce9b5ae46b8675049b5bc4d9525a0 \ No newline at end of file From c7c5e671cfdca54cb9ed93c9f8c19da7aebc7636 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Oct 2016 16:26:46 -0400 Subject: [PATCH 15/32] Use ordered dict in test_multipole --- tests/test_multipole/inputs_true.dat | 2 +- tests/test_multipole/test_multipole.py | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 3b2b43ec7..83e9fc6fd 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -1 +1 @@ -2be927608035759f52a2ac88b2c84e28c06e0f7905187bfd4695fecd80f417e727d8879c52fe03253938f7aa0301061f72dce9b5ae46b8675049b5bc4d9525a0 \ No newline at end of file +28e3b5cba12061ff7c410585117d40425d830a53fa2e7324d8f96deee28a57fb4cc4de393aa600f44816e72a0f66f22f5f4c65d52977d99f6a0d53275d263067 \ No newline at end of file diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index 32b38dd76..035a6c894 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from collections import OrderedDict import os import sys sys.path.insert(0, os.pardir) @@ -68,7 +69,9 @@ class MultipoleTestHarness(PyAPITestHarness): sets_file.particles = 1000 sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1])) sets_file.output = {'summary': True} - sets_file.temperature = {'tolerance': 1000, 'multipole': True} + sets_file.temperature = OrderedDict() + sets_file.temperature['tolerance'] = 1000 + sets_file.temperature['multipole'] = True sets_file.export_to_xml() #################### From e96849c3842142641053ba3619ff202fbd28ead5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Oct 2016 17:02:49 -0400 Subject: [PATCH 16/32] Improve 0K substitution warning messages --- src/input_xml.F90 | 4 ++-- src/nuclide_header.F90 | 24 +++++++++++++++--------- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 16892397a..5aa22fba5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -5765,7 +5765,7 @@ contains file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) call nuclides(i_nuclide) % from_hdf5(group_id, nuc_temps(i_nuclide), & - temperature_method, temperature_tolerance) + temperature_method, temperature_tolerance, master) call close_group(group_id) call file_close(file_id) @@ -5999,7 +5999,7 @@ contains group_id = open_group(file_id, name) method = TEMPERATURE_NEAREST call resonant_nuc % from_hdf5(group_id, temperature, & - method, 1000.0_8) + method, 1000.0_8, master) call close_group(group_id) call file_close(file_id) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 4235f4234..93832869c 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -186,12 +186,16 @@ module nuclide_header end subroutine nuclide_clear - subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance) - class(Nuclide), intent(inout) :: this - integer(HID_T), intent(in) :: group_id + subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance, & + master) + class(Nuclide), intent(inout) :: this + integer(HID_T), intent(in) :: group_id type(VectorReal), intent(in) :: temperature ! list of desired temperatures - integer, intent(inout) :: method - real(8), intent(in) :: tolerance + integer, intent(inout) :: method + real(8), intent(in) :: tolerance + logical, intent(in) :: master ! if this is the master proc + ! (can't use global cuz + ! circular dependance) integer :: i integer :: storage_type @@ -265,15 +269,17 @@ module nuclide_header call temps_to_read % push_back(nint(temp_actual)) ! Write warning for resonance scattering data if 0K is not available - if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0) then + if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0 & + .and. master) then call warning(trim(this % name) // " does not contain 0K data & &needed for resonance scattering options selected. Using & - &data at " // trim(to_str(nint(temp_actual))) // " K instead.") + &data at " // trim(to_str(temp_actual)) & + // " K instead.") end if end if else - call fatal_error("Nuclear data library does not contain cross sections & - &for " // trim(this % name) // " at or near " // & + call fatal_error("Nuclear data library does not contain cross & + §ions for " // trim(this % name) // " at or near " // & trim(to_str(nint(temp_desired))) // " K.") end if end do From b4df08ac336c3a26c5a94a55ecba2ef006e7cb36 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Oct 2016 20:54:57 -0400 Subject: [PATCH 17/32] Sort the output of temperature dictionary --- openmc/settings.py | 2 +- tests/test_multipole/inputs_true.dat | 2 +- tests/test_multipole/test_multipole.py | 5 +---- 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 84fe20485..335de839a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1053,7 +1053,7 @@ class Settings(object): def _create_temperature_subelements(self): if self.temperature: - for key, value in self.temperature.items(): + for key, value in sorted(self.temperature.items()): element = ET.SubElement(self._settings_file, "temperature_{}".format(key)) element.text = str(value) diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 83e9fc6fd..3b2b43ec7 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -1 +1 @@ -28e3b5cba12061ff7c410585117d40425d830a53fa2e7324d8f96deee28a57fb4cc4de393aa600f44816e72a0f66f22f5f4c65d52977d99f6a0d53275d263067 \ No newline at end of file +2be927608035759f52a2ac88b2c84e28c06e0f7905187bfd4695fecd80f417e727d8879c52fe03253938f7aa0301061f72dce9b5ae46b8675049b5bc4d9525a0 \ No newline at end of file diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index 035a6c894..32b38dd76 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -1,5 +1,4 @@ #!/usr/bin/env python -from collections import OrderedDict import os import sys sys.path.insert(0, os.pardir) @@ -69,9 +68,7 @@ class MultipoleTestHarness(PyAPITestHarness): sets_file.particles = 1000 sets_file.source = Source(space=Box([-1, -1, -1], [1, 1, 1])) sets_file.output = {'summary': True} - sets_file.temperature = OrderedDict() - sets_file.temperature['tolerance'] = 1000 - sets_file.temperature['multipole'] = True + sets_file.temperature = {'tolerance': 1000, 'multipole': True} sets_file.export_to_xml() #################### From 9a8d42baf32dfa19df0f704cc1fe23f4a719ec8e Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 8 Oct 2016 10:40:28 -0400 Subject: [PATCH 18/32] modified travis.yml --- .travis.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 46bb847e6..ed81e6716 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,4 +1,5 @@ -sudo: false +sudo: required +dist: trusty language: python python: - "2.7" From 73d9f7ce019f700211c92cf1ca4ae1abd7b2ec35 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 8 Oct 2016 11:20:13 -0400 Subject: [PATCH 19/32] removed low-hanging gnu fortran 4.6 specific code --- src/global.F90 | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index bea5f61a8..d971581ba 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -468,12 +468,7 @@ contains do i = 1, size(nuclides) call nuclides(i) % clear() end do - - ! WARNING: The following statement should work but doesn't under gfortran - ! 4.6 because of a bug. Technically, commenting this out leaves a memory - ! leak. - - ! deallocate(nuclides) + deallocate(nuclides) end if if (allocated(nuclides_0K)) then From 490a0d8338dc2b1e24ed2f39961a121618ce1e79 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 8 Oct 2016 11:34:28 -0400 Subject: [PATCH 20/32] Updated reference to gfrotran 4.6.0 in the docs as being the minimum version --- docs/source/usersguide/install.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 3dab08e9d..925008495 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -51,7 +51,7 @@ Prerequisites installed on your machine. Since a number of Fortran 2003/2008 features are used in the code, it is recommended that you use the latest version of whatever compiler you choose. For gfortran_, it is necessary to use - version 4.6.0 or above. + version 4.8.0 or above. If you are using Debian or a Debian derivative such as Ubuntu, you can install the gfortran compiler using the following command:: From 7eefb7306ff0f8745af1447789d8cddc9778617f Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 9 Oct 2016 15:10:24 -0400 Subject: [PATCH 21/32] Fixed errors in tally.F90 which would have resulted in erroneous tallied values. --- src/tally.F90 | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 3c008e30e..f7d8ab5a2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -149,7 +149,7 @@ contains if (survival_biasing) then ! We need to account for the fact that some weight was already ! absorbed - score = p % last_wgt + p % absorb_wgt * flux + score = (p % last_wgt + p % absorb_wgt) * flux else score = p % last_wgt * flux end if @@ -417,6 +417,13 @@ contains case (SCORE_PROMPT_NU_FISSION) + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if + if (t % estimator == ESTIMATOR_ANALOG) then if (survival_biasing .or. p % fission) then if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -453,13 +460,6 @@ contains end if else - ! make sure the correct energy is used - if (t % estimator == ESTIMATOR_TRACKLENGTH) then - E = p % E - else - E = p % last_E - end if - if (i_nuclide > 0) then score = micro_xs(i_nuclide) % fission * nuclides(i_nuclide) % & nu(E, EMISSION_PROMPT) * atom_density * flux @@ -675,6 +675,12 @@ contains case (SCORE_DECAY_RATE) + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if ! Set the delayedgroup filter index dg_filter = t % find_filter(FILTER_DELAYEDGROUP) From 65513b2664b4686f15c39592b77d6f4d94098ed3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 10 Oct 2016 09:04:57 -0400 Subject: [PATCH 22/32] Changes to testing scripts to less over-subscribe Travis (3 OMP threads to 2 and 3 MPI processes to 2) --- .travis.yml | 2 +- tests/testing_harness.py | 2 +- tests/travis_install.sh | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index ed81e6716..4236c36d1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -55,6 +55,6 @@ before_script: script: - cd tests - - export OMP_NUM_THREADS=3 + - export OMP_NUM_THREADS=2 - ./travis.sh - cd .. diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 530013164..150d124c8 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -24,7 +24,7 @@ class TestHarness(object): self.parser = OptionParser() self.parser.add_option('--exe', dest='exe', default='openmc') self.parser.add_option('--mpi_exec', dest='mpi_exec', default=None) - self.parser.add_option('--mpi_np', dest='mpi_np', type=int, default=3) + self.parser.add_option('--mpi_np', dest='mpi_np', type=int, default=2) self.parser.add_option('--update', dest='update', action='store_true', default=False) self._opts = None diff --git a/tests/travis_install.sh b/tests/travis_install.sh index df1c81314..2e35c4f65 100755 --- a/tests/travis_install.sh +++ b/tests/travis_install.sh @@ -8,7 +8,7 @@ if [[ ! -e $HOME/mpich_install/bin/mpiexec ]]; then tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1 cd mpich-3.1.3 ./configure --prefix=$HOME/mpich_install -q - make >/dev/null 2>&1 + make -j 2 >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. fi @@ -22,7 +22,7 @@ if [[ ! -e $HOME/phdf5_install/bin/h5pfc ]]; then ./configure \ --prefix=$HOME/phdf5_install -q --enable-fortran \ --enable-fortran2003 --enable-parallel - make >/dev/null 2>&1 + make -j 2 >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. fi @@ -33,7 +33,7 @@ if [[ ! -e $HOME/hdf5_install/bin/h5fc ]]; then cd hdf5-1.8.15 CC=gcc FC=gfortran ./configure --prefix=$HOME/hdf5_install -q \ --enable-fortran --enable-fortran2003 - make -j >/dev/null 2>&1 + make -j 2 >/dev/null 2>&1 make install >/dev/null 2>&1 cd .. fi From 4f3454edead1eacb3aed4b53e3aa9ff79724b316 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 10 Oct 2016 09:49:01 -0400 Subject: [PATCH 23/32] Significantly reduced the number of nuclides in the standard (HM problem I believe) input set, reducing the total test time by a good amount --- tests/input_set.py | 96 ----- tests/test_asymmetric_lattice/inputs_true.dat | 2 +- .../test_asymmetric_lattice/results_true.dat | 2 +- tests/test_filter_mesh/inputs_true.dat | 2 +- tests/test_filter_mesh/results_true.dat | 2 +- tests/test_iso_in_lab/inputs_true.dat | 2 +- tests/test_iso_in_lab/results_true.dat | 2 +- tests/test_mgxs_library_mesh/inputs_true.dat | 2 +- tests/test_mgxs_library_mesh/results_true.dat | 348 +++++++++--------- tests/test_tallies/inputs_true.dat | 2 +- tests/test_tallies/results_true.dat | 2 +- tests/test_tally_aggregation/inputs_true.dat | 2 +- tests/test_tally_aggregation/results_true.dat | 2 +- .../test_tally_aggregation.py | 6 +- tests/test_tally_arithmetic/inputs_true.dat | 2 +- .../test_tally_arithmetic.py | 4 +- tests/test_tally_slice_merge/inputs_true.dat | 2 +- tests/test_tally_slice_merge/results_true.dat | 96 ++--- 18 files changed, 240 insertions(+), 336 deletions(-) diff --git a/tests/input_set.py b/tests/input_set.py index 94576acf0..f3f9e88bd 100644 --- a/tests/input_set.py +++ b/tests/input_set.py @@ -27,37 +27,8 @@ class InputSet(object): fuel.set_density('g/cm3', 10.062) fuel.add_nuclide("U234", 4.9476e-6) fuel.add_nuclide("U235", 4.8218e-4) - fuel.add_nuclide("U236", 9.0402e-5) fuel.add_nuclide("U238", 2.1504e-2) - fuel.add_nuclide("Np237", 7.3733e-6) - fuel.add_nuclide("Pu238", 1.5148e-6) - fuel.add_nuclide("Pu239", 1.3955e-4) - fuel.add_nuclide("Pu240", 3.4405e-5) - fuel.add_nuclide("Pu241", 2.1439e-5) - fuel.add_nuclide("Pu242", 3.7422e-6) - fuel.add_nuclide("Am241", 4.5041e-7) - fuel.add_nuclide("Am242_m1", 9.2301e-9) - fuel.add_nuclide("Am243", 4.7878e-7) - fuel.add_nuclide("Cm242", 1.0485e-7) - fuel.add_nuclide("Cm243", 1.4268e-9) - fuel.add_nuclide("Cm244", 8.8756e-8) - fuel.add_nuclide("Cm245", 3.5285e-9) - fuel.add_nuclide("Mo95", 2.6497e-5) - fuel.add_nuclide("Tc99", 3.2772e-5) - fuel.add_nuclide("Ru101", 3.0742e-5) - fuel.add_nuclide("Ru103", 2.3505e-6) - fuel.add_nuclide("Ag109", 2.0009e-6) fuel.add_nuclide("Xe135", 1.0801e-8) - fuel.add_nuclide("Cs133", 3.4612e-5) - fuel.add_nuclide("Nd143", 2.6078e-5) - fuel.add_nuclide("Nd145", 1.9898e-5) - fuel.add_nuclide("Sm147", 1.6128e-6) - fuel.add_nuclide("Sm149", 1.1627e-7) - fuel.add_nuclide("Sm150", 7.1727e-6) - fuel.add_nuclide("Sm151", 5.4947e-7) - fuel.add_nuclide("Sm152", 3.0221e-6) - fuel.add_nuclide("Eu153", 2.6209e-6) - fuel.add_nuclide("Gd155", 1.5369e-9) fuel.add_nuclide("O16", 4.5737e-2) clad = openmc.Material(name='Cladding', material_id=2) @@ -93,27 +64,10 @@ class InputSet(object): rpv_steel.add_nuclide("Fe58", 0.00282159, 'wo') rpv_steel.add_nuclide("Ni58", 0.0067198, 'wo') rpv_steel.add_nuclide("Ni60", 0.0026776, 'wo') - rpv_steel.add_nuclide("Ni61", 0.0001183, 'wo') - rpv_steel.add_nuclide("Ni62", 0.0003835, 'wo') - rpv_steel.add_nuclide("Ni64", 0.0001008, 'wo') rpv_steel.add_nuclide("Mn55", 0.01, 'wo') - rpv_steel.add_nuclide("Mo92", 0.000849, 'wo') - rpv_steel.add_nuclide("Mo94", 0.0005418, 'wo') - rpv_steel.add_nuclide("Mo95", 0.0009438, 'wo') - rpv_steel.add_nuclide("Mo96", 0.0010002, 'wo') - rpv_steel.add_nuclide("Mo97", 0.0005796, 'wo') - rpv_steel.add_nuclide("Mo98", 0.0014814, 'wo') - rpv_steel.add_nuclide("Mo100", 0.0006042, 'wo') - rpv_steel.add_nuclide("Si28", 0.00367464, 'wo') - rpv_steel.add_nuclide("Si29", 0.00019336, 'wo') - rpv_steel.add_nuclide("Si30", 0.000132, 'wo') - rpv_steel.add_nuclide("Cr50", 0.00010435, 'wo') rpv_steel.add_nuclide("Cr52", 0.002092475, 'wo') - rpv_steel.add_nuclide("Cr53", 0.00024185, 'wo') - rpv_steel.add_nuclide("Cr54", 6.1325e-05, 'wo') rpv_steel.add_nuclide("C0", 0.0025, 'wo') rpv_steel.add_nuclide("Cu63", 0.0013696, 'wo') - rpv_steel.add_nuclide("Cu65", 0.0006304, 'wo') lower_rad_ref = openmc.Material(name='Lower radial reflector', material_id=6) @@ -127,18 +81,8 @@ class InputSet(object): lower_rad_ref.add_nuclide("Fe57", 0.01362750048, 'wo') lower_rad_ref.add_nuclide("Fe58", 0.001848545204, 'wo') lower_rad_ref.add_nuclide("Ni58", 0.055298376566, 'wo') - lower_rad_ref.add_nuclide("Ni60", 0.022034425592, 'wo') - lower_rad_ref.add_nuclide("Ni61", 0.000973510811, 'wo') - lower_rad_ref.add_nuclide("Ni62", 0.003155886695, 'wo') - lower_rad_ref.add_nuclide("Ni64", 0.000829500336, 'wo') lower_rad_ref.add_nuclide("Mn55", 0.0182870, 'wo') - lower_rad_ref.add_nuclide("Si28", 0.00839976771, 'wo') - lower_rad_ref.add_nuclide("Si29", 0.00044199679, 'wo') - lower_rad_ref.add_nuclide("Si30", 0.0003017355, 'wo') - lower_rad_ref.add_nuclide("Cr50", 0.007251360806, 'wo') lower_rad_ref.add_nuclide("Cr52", 0.145407678031, 'wo') - lower_rad_ref.add_nuclide("Cr53", 0.016806340306, 'wo') - lower_rad_ref.add_nuclide("Cr54", 0.004261520857, 'wo') lower_rad_ref.add_s_alpha_beta('c_H_in_H2O') upper_rad_ref = openmc.Material(name='Upper radial reflector /' @@ -153,18 +97,8 @@ class InputSet(object): upper_rad_ref.add_nuclide("Fe57", 0.01375486056, 'wo') upper_rad_ref.add_nuclide("Fe58", 0.001865821363, 'wo') upper_rad_ref.add_nuclide("Ni58", 0.055815129186, 'wo') - upper_rad_ref.add_nuclide("Ni60", 0.022240333032, 'wo') - upper_rad_ref.add_nuclide("Ni61", 0.000982608081, 'wo') - upper_rad_ref.add_nuclide("Ni62", 0.003185377845, 'wo') - upper_rad_ref.add_nuclide("Ni64", 0.000837251856, 'wo') upper_rad_ref.add_nuclide("Mn55", 0.0184579, 'wo') - upper_rad_ref.add_nuclide("Si28", 0.00847831314, 'wo') - upper_rad_ref.add_nuclide("Si29", 0.00044612986, 'wo') - upper_rad_ref.add_nuclide("Si30", 0.000304557, 'wo') - upper_rad_ref.add_nuclide("Cr50", 0.00731912987, 'wo') upper_rad_ref.add_nuclide("Cr52", 0.146766614995, 'wo') - upper_rad_ref.add_nuclide("Cr53", 0.01696340737, 'wo') - upper_rad_ref.add_nuclide("Cr54", 0.004301347765, 'wo') upper_rad_ref.add_s_alpha_beta('c_H_in_H2O') bot_plate = openmc.Material(name='Bottom plate region', material_id=8) @@ -178,18 +112,8 @@ class InputSet(object): bot_plate.add_nuclide("Fe57", 0.014750478, 'wo') bot_plate.add_nuclide("Fe58", 0.002000875025, 'wo') bot_plate.add_nuclide("Ni58", 0.059855207342, 'wo') - bot_plate.add_nuclide("Ni60", 0.023850159704, 'wo') - bot_plate.add_nuclide("Ni61", 0.001053732407, 'wo') - bot_plate.add_nuclide("Ni62", 0.003415945715, 'wo') - bot_plate.add_nuclide("Ni64", 0.000897854832, 'wo') bot_plate.add_nuclide("Mn55", 0.0197940, 'wo') - bot_plate.add_nuclide("Si28", 0.00909197802, 'wo') - bot_plate.add_nuclide("Si29", 0.00047842098, 'wo') - bot_plate.add_nuclide("Si30", 0.000326601, 'wo') - bot_plate.add_nuclide("Cr50", 0.007848910646, 'wo') bot_plate.add_nuclide("Cr52", 0.157390026871, 'wo') - bot_plate.add_nuclide("Cr53", 0.018191270146, 'wo') - bot_plate.add_nuclide("Cr54", 0.004612692337, 'wo') bot_plate.add_s_alpha_beta('c_H_in_H2O') bot_nozzle = openmc.Material(name='Bottom nozzle region', @@ -204,18 +128,8 @@ class InputSet(object): bot_nozzle.add_nuclide("Fe57", 0.01163454624, 'wo') bot_nozzle.add_nuclide("Fe58", 0.001578204652, 'wo') bot_nozzle.add_nuclide("Ni58", 0.047211231662, 'wo') - bot_nozzle.add_nuclide("Ni60", 0.018811987544, 'wo') - bot_nozzle.add_nuclide("Ni61", 0.000831139127, 'wo') - bot_nozzle.add_nuclide("Ni62", 0.002694352115, 'wo') - bot_nozzle.add_nuclide("Ni64", 0.000708189552, 'wo') bot_nozzle.add_nuclide("Mn55", 0.0156126, 'wo') - bot_nozzle.add_nuclide("Si28", 0.007171335558, 'wo') - bot_nozzle.add_nuclide("Si29", 0.000377356542, 'wo') - bot_nozzle.add_nuclide("Si30", 0.0002576079, 'wo') - bot_nozzle.add_nuclide("Cr50", 0.006190885148, 'wo') bot_nozzle.add_nuclide("Cr52", 0.124142524198, 'wo') - bot_nozzle.add_nuclide("Cr53", 0.014348496148, 'wo') - bot_nozzle.add_nuclide("Cr54", 0.003638294506, 'wo') bot_nozzle.add_s_alpha_beta('c_H_in_H2O') top_nozzle = openmc.Material(name='Top nozzle region', material_id=10) @@ -229,18 +143,8 @@ class InputSet(object): top_nozzle.add_nuclide("Fe57", 0.0101152584, 'wo') top_nozzle.add_nuclide("Fe58", 0.00137211607, 'wo') top_nozzle.add_nuclide("Ni58", 0.04104621835, 'wo') - top_nozzle.add_nuclide("Ni60", 0.0163554502, 'wo') - top_nozzle.add_nuclide("Ni61", 0.000722605975, 'wo') - top_nozzle.add_nuclide("Ni62", 0.002342513875, 'wo') - top_nozzle.add_nuclide("Ni64", 0.0006157116, 'wo') top_nozzle.add_nuclide("Mn55", 0.0135739, 'wo') - top_nozzle.add_nuclide("Si28", 0.006234853554, 'wo') - top_nozzle.add_nuclide("Si29", 0.000328078746, 'wo') - top_nozzle.add_nuclide("Si30", 0.0002239677, 'wo') - top_nozzle.add_nuclide("Cr50", 0.005382452306, 'wo') top_nozzle.add_nuclide("Cr52", 0.107931450781, 'wo') - top_nozzle.add_nuclide("Cr53", 0.012474806806, 'wo') - top_nozzle.add_nuclide("Cr54", 0.003163190107, 'wo') top_nozzle.add_s_alpha_beta('c_H_in_H2O') top_fa = openmc.Material(name='Top of fuel assemblies', material_id=11) diff --git a/tests/test_asymmetric_lattice/inputs_true.dat b/tests/test_asymmetric_lattice/inputs_true.dat index 9d278724f..7e48dcf43 100644 --- a/tests/test_asymmetric_lattice/inputs_true.dat +++ b/tests/test_asymmetric_lattice/inputs_true.dat @@ -1 +1 @@ -dfb59bace10a91bb7ffc871d8ee87e91d94754bb8bb002ac6088f80fe0f480741c0489f74b753fc37158d0ff0f1368739ea60638b42083791311eefeac79168e \ No newline at end of file +6bcc9cca24d42995bdff9bf9aca5e852c2dbca5cfb42a12ac637def9cf5cac227654182fc9cf9e17d07cf2e9af11fea832e3ae0eb7001cc09856f73d219664f9 \ No newline at end of file diff --git a/tests/test_asymmetric_lattice/results_true.dat b/tests/test_asymmetric_lattice/results_true.dat index a33b9c9e5..c753afbc7 100644 --- a/tests/test_asymmetric_lattice/results_true.dat +++ b/tests/test_asymmetric_lattice/results_true.dat @@ -1 +1 @@ -bc8bef8121f9b6470e4fea817a4e48eabb1ecba1f42761a4cbd77d71181bf9e1612df4a3d6ddfbcd08a3086ac873e5f3c3e560bf96b2b7c959a2f7aad7e4e08d \ No newline at end of file +a2848fdb0a12c99ce31f4ddee766e0cf33bd5c6feca927bcb4bceca6fbb094bb1309fb8548589a99fcb09a830911457b9175b460aa45773822f8112a34b3b5c1 \ No newline at end of file diff --git a/tests/test_filter_mesh/inputs_true.dat b/tests/test_filter_mesh/inputs_true.dat index b683ebaeb..d98e2f6f7 100644 --- a/tests/test_filter_mesh/inputs_true.dat +++ b/tests/test_filter_mesh/inputs_true.dat @@ -1 +1 @@ -5a9e65b8a8c9d7c575fc48c5d289bbc805739729a33d83aa79985473d83c8cc3a0c7dad8e95221090917c35ac4842667c8c9daecd06dc6b909a92475f5083753 \ No newline at end of file +c0882d16048d434219d32ed7b875615e00142b55ab78be5b51e0195e23cb534a3aa5a8e9c212f0bfff96c452177315deb1264421d688438814560b4131c88930 \ No newline at end of file diff --git a/tests/test_filter_mesh/results_true.dat b/tests/test_filter_mesh/results_true.dat index 2cc3c9532..06cd450ef 100644 --- a/tests/test_filter_mesh/results_true.dat +++ b/tests/test_filter_mesh/results_true.dat @@ -1 +1 @@ -89387dd9e5b962c773e1782ff829846968dc456d899963f5dd98761fbde73a3f81676717ff3599bc26a1b77247826c38756735a9b2b329b80ad66286a62bd3f9 \ No newline at end of file +041b7e79d384771bad5bc6800138bd20be6af7d6e1543bdcfef521de3df33ef7b7aca00ea3ccb0234c3a4d7d738dfe17987b5c1a59474efef630217f8ee1af1e \ No newline at end of file diff --git a/tests/test_iso_in_lab/inputs_true.dat b/tests/test_iso_in_lab/inputs_true.dat index 310bccb13..454f891cd 100644 --- a/tests/test_iso_in_lab/inputs_true.dat +++ b/tests/test_iso_in_lab/inputs_true.dat @@ -1 +1 @@ -4b3d0270a479e65579b305d1c2339b76971790bc7371c685efa6e2d341980fec301cf0859c34796ae04ae98aa01ab8b4905a8d8a3a916895c36d82ca6b58fb39 \ No newline at end of file +57d1ece4aa9633e5fc6d2fce9f7bba444d02046e0cdadc6e53efdf591b6e50fd0e5a14b956b7ad10710b4535b10354899e2ffc783fd78de5ccd09a061cb5678d \ No newline at end of file diff --git a/tests/test_iso_in_lab/results_true.dat b/tests/test_iso_in_lab/results_true.dat index 354ccb0f8..5b2c5dd00 100644 --- a/tests/test_iso_in_lab/results_true.dat +++ b/tests/test_iso_in_lab/results_true.dat @@ -1,2 +1,2 @@ k-combined: -9.638450E-01 1.237705E-02 +9.753410E-01 6.608245E-02 diff --git a/tests/test_mgxs_library_mesh/inputs_true.dat b/tests/test_mgxs_library_mesh/inputs_true.dat index 6e3fc6e63..a87549a60 100644 --- a/tests/test_mgxs_library_mesh/inputs_true.dat +++ b/tests/test_mgxs_library_mesh/inputs_true.dat @@ -1 +1 @@ -88b1a1b52bcae466bab63038d03e9b2fb9f2a28aac3d73573ca00d5f0c116aa9e1c97341215b80266e8317a0549f11d9adcb8726c82fcc1741a5a730e10c8186 \ No newline at end of file +5011f0aec7bb04e22f96cec2c23ae485fe21371dc03e16cb38695b3765cb92545459a897184b501bb99f8e7d162cfbf85eb3fb52ce2d82b2d2ca722b567cbff3 \ No newline at end of file diff --git a/tests/test_mgxs_library_mesh/results_true.dat b/tests/test_mgxs_library_mesh/results_true.dat index 8ad0f2eee..16d20391c 100644 --- a/tests/test_mgxs_library_mesh/results_true.dat +++ b/tests/test_mgxs_library_mesh/results_true.dat @@ -1,168 +1,168 @@ mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.640786 0.044177 -1 1 2 1 1 total 0.615276 0.104046 -2 2 1 1 1 total 0.660597 0.128423 -3 2 2 1 1 total 0.646999 0.186709 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.36665 0.048814 -1 1 2 1 1 total 0.36356 0.074111 -2 2 1 1 1 total 0.40784 0.096486 -3 2 2 1 1 total 0.41456 0.160443 +0 1 1 1 1 total 0.654966 0.098415 +1 1 2 1 1 total 0.639357 0.282107 +2 2 1 1 1 total 0.713534 0.079789 +3 2 2 1 1 total 0.641095 0.091519 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.366650 0.048814 -1 1 2 1 1 total 0.363560 0.074111 -2 2 1 1 1 total 0.407840 0.096486 -3 2 2 1 1 total 0.414593 0.160436 +0 1 1 1 1 total 0.413423 0.087250 +1 1 2 1 1 total 0.392074 0.244272 +2 2 1 1 1 total 0.458841 0.087921 +3 2 2 1 1 total 0.403898 0.074343 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.025749 0.002863 -1 1 2 1 1 total 0.022988 0.004099 -2 2 1 1 1 total 0.028400 0.005275 -3 2 2 1 1 total 0.027589 0.010350 +0 1 1 1 1 total 0.413423 0.087250 +1 1 2 1 1 total 0.392074 0.244272 +2 2 1 1 1 total 0.458841 0.087921 +3 2 2 1 1 total 0.403898 0.074343 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.015861 0.002876 -1 1 2 1 1 total 0.014403 0.003542 -2 2 1 1 1 total 0.017280 0.004371 -3 2 2 1 1 total 0.018061 0.010110 +0 1 1 1 1 total 0.021476 0.004248 +1 1 2 1 1 total 0.020653 0.008355 +2 2 1 1 1 total 0.027384 0.003568 +3 2 2 1 1 total 0.021826 0.004584 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.009888 0.001077 -1 1 2 1 1 total 0.008585 0.001552 -2 2 1 1 1 total 0.011121 0.002456 -3 2 2 1 1 total 0.009527 0.003659 +0 1 1 1 1 total 0.013234 0.004146 +1 1 2 1 1 total 0.012342 0.006800 +2 2 1 1 1 total 0.016807 0.003428 +3 2 2 1 1 total 0.013253 0.004795 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.026065 0.002907 -1 1 2 1 1 total 0.022596 0.004062 -2 2 1 1 1 total 0.029084 0.006430 -3 2 2 1 1 total 0.025066 0.009687 +0 1 1 1 1 total 0.008241 0.001798 +1 1 2 1 1 total 0.008311 0.003296 +2 2 1 1 1 total 0.010577 0.001333 +3 2 2 1 1 total 0.008573 0.002017 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 1.938476 0.211550 -1 1 2 1 1 total 1.682799 0.303764 -2 2 1 1 1 total 2.177360 0.480780 -3 2 2 1 1 total 1.864890 0.715661 +0 1 1 1 1 total 0.020322 0.004415 +1 1 2 1 1 total 0.020546 0.008145 +2 2 1 1 1 total 0.026008 0.003213 +3 2 2 1 1 total 0.021015 0.004911 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.615037 0.041754 -1 1 2 1 1 total 0.592288 0.100439 -2 2 1 1 1 total 0.632196 0.123878 -3 2 2 1 1 total 0.619410 0.177190 +0 1 1 1 1 total 1.596880 0.348354 +1 1 2 1 1 total 1.610266 0.638515 +2 2 1 1 1 total 2.048209 0.257682 +3 2 2 1 1 total 1.660283 0.390011 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.584014 0.054315 -1 1 2 1 1 total 0.587256 0.084833 -2 2 1 1 1 total 0.622514 0.111323 -3 2 2 1 1 total 0.613792 0.168612 +0 1 1 1 1 total 0.633490 0.094536 +1 1 2 1 1 total 0.618705 0.273934 +2 2 1 1 1 total 0.686150 0.076703 +3 2 2 1 1 total 0.619269 0.087616 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.638348 0.097266 +1 1 2 1 1 total 0.588378 0.266032 +2 2 1 1 1 total 0.698373 0.092394 +3 2 2 1 1 total 0.625643 0.075089 mesh 1 group in group out nuclide moment mean std. dev. x y z -0 1 1 1 1 1 total P0 0.584014 0.054315 -1 1 1 1 1 1 total P1 0.243427 0.025488 -2 1 1 1 1 1 total P2 0.089236 0.007357 -3 1 1 1 1 1 total P3 0.008994 0.005768 -4 1 2 1 1 1 total P0 0.587256 0.084833 -5 1 2 1 1 1 total P1 0.245120 0.041033 -6 1 2 1 1 1 total P2 0.086784 0.016255 -7 1 2 1 1 1 total P3 0.008660 0.004755 -8 2 1 1 1 1 total P0 0.622514 0.111323 -9 2 1 1 1 1 total P1 0.239376 0.042594 -10 2 1 1 1 1 total P2 0.088386 0.017200 -11 2 1 1 1 1 total P3 -0.001243 0.005639 -12 2 2 1 1 1 total P0 0.612950 0.167940 -13 2 2 1 1 1 total P1 0.226176 0.061882 -14 2 2 1 1 1 total P2 0.086593 0.026126 -15 2 2 1 1 1 total P3 0.009672 0.011995 +0 1 1 1 1 1 total P0 0.638348 0.097266 +1 1 1 1 1 1 total P1 0.247099 0.035574 +2 1 1 1 1 1 total P2 0.092195 0.010891 +3 1 1 1 1 1 total P3 0.013224 0.004528 +4 1 2 1 1 1 total P0 0.588378 0.266032 +5 1 2 1 1 1 total P1 0.221552 0.102564 +6 1 2 1 1 1 total P2 0.069798 0.034699 +7 1 2 1 1 1 total P3 -0.003039 0.009633 +8 2 1 1 1 1 total P0 0.698373 0.092394 +9 2 1 1 1 1 total P1 0.261835 0.035253 +10 2 1 1 1 1 total P2 0.096206 0.013947 +11 2 1 1 1 1 total P3 0.016973 0.004808 +12 2 2 1 1 1 total P0 0.625643 0.075089 +13 2 2 1 1 1 total P1 0.244646 0.028604 +14 2 2 1 1 1 total P2 0.088580 0.012369 +15 2 2 1 1 1 total P3 0.019989 0.013950 mesh 1 group in group out nuclide moment mean std. dev. x y z -0 1 1 1 1 1 total P0 0.584014 0.054315 -1 1 1 1 1 1 total P1 0.243427 0.025488 -2 1 1 1 1 1 total P2 0.089236 0.007357 -3 1 1 1 1 1 total P3 0.008994 0.005768 -4 1 2 1 1 1 total P0 0.587256 0.084833 -5 1 2 1 1 1 total P1 0.245120 0.041033 -6 1 2 1 1 1 total P2 0.086784 0.016255 -7 1 2 1 1 1 total P3 0.008660 0.004755 -8 2 1 1 1 1 total P0 0.622514 0.111323 -9 2 1 1 1 1 total P1 0.239376 0.042594 -10 2 1 1 1 1 total P2 0.088386 0.017200 -11 2 1 1 1 1 total P3 -0.001243 0.005639 -12 2 2 1 1 1 total P0 0.613792 0.168612 -13 2 2 1 1 1 total P1 0.226142 0.061856 -14 2 2 1 1 1 total P2 0.086174 0.025979 -15 2 2 1 1 1 total P3 0.009721 0.012027 +0 1 1 1 1 1 total P0 0.638348 0.097266 +1 1 1 1 1 1 total P1 0.247099 0.035574 +2 1 1 1 1 1 total P2 0.092195 0.010891 +3 1 1 1 1 1 total P3 0.013224 0.004528 +4 1 2 1 1 1 total P0 0.588378 0.266032 +5 1 2 1 1 1 total P1 0.221552 0.102564 +6 1 2 1 1 1 total P2 0.069798 0.034699 +7 1 2 1 1 1 total P3 -0.003039 0.009633 +8 2 1 1 1 1 total P0 0.698373 0.092394 +9 2 1 1 1 1 total P1 0.261835 0.035253 +10 2 1 1 1 1 total P2 0.096206 0.013947 +11 2 1 1 1 1 total P3 0.016973 0.004808 +12 2 2 1 1 1 total P0 0.625643 0.075089 +13 2 2 1 1 1 total P1 0.244646 0.028604 +14 2 2 1 1 1 total P2 0.088580 0.012369 +15 2 2 1 1 1 total P3 0.019989 0.013950 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 1.0 0.153265 +1 1 2 1 1 1 total 1.0 0.454973 +2 2 1 1 1 1 total 1.0 0.146747 +3 2 2 1 1 1 total 1.0 0.141824 mesh 1 group in group out nuclide mean std. dev. x y z -0 1 1 1 1 1 total 1.000000 0.088094 -1 1 2 1 1 1 total 1.000000 0.126864 -2 2 1 1 1 1 total 1.000000 0.160891 -3 2 2 1 1 1 total 1.001374 0.305883 - mesh 1 group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.027395 0.004680 -1 1 2 1 1 1 total 0.019384 0.002846 -2 2 1 1 1 1 total 0.022914 0.006025 -3 2 2 1 1 1 total 0.029629 0.006292 +0 1 1 1 1 1 total 0.021059 0.003031 +1 1 2 1 1 1 total 0.017348 0.008786 +2 2 1 1 1 1 total 0.020409 0.003354 +3 2 2 1 1 1 total 0.011105 0.003806 mesh 1 group out nuclide mean std. dev. x y z -0 1 1 1 1 total 1.0 0.220956 -1 1 2 1 1 total 1.0 0.132140 -2 2 1 1 1 total 1.0 0.316565 -3 2 2 1 1 total 1.0 0.181577 +0 1 1 1 1 total 1.0 0.135958 +1 1 2 1 1 total 1.0 0.557756 +2 2 1 1 1 total 1.0 0.201340 +3 2 2 1 1 total 1.0 0.475608 mesh 1 group out nuclide mean std. dev. x y z -0 1 1 1 1 total 1.0 0.222246 -1 1 2 1 1 total 1.0 0.132140 -2 2 1 1 1 total 1.0 0.316565 -3 2 2 1 1 total 1.0 0.181577 +0 1 1 1 1 total 1.0 0.133151 +1 1 2 1 1 total 1.0 0.557756 +2 2 1 1 1 total 1.0 0.201340 +3 2 2 1 1 total 1.0 0.475608 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 3.610522e-07 3.169931e-08 -1 1 2 1 1 total 3.097784e-07 5.252025e-08 -2 2 1 1 1 total 3.942353e-07 8.459167e-08 -3 2 2 1 1 total 3.799163e-07 1.806470e-07 +0 1 1 1 1 total 4.697405e-07 8.304381e-08 +1 1 2 1 1 total 4.173069e-07 1.694647e-07 +2 2 1 1 1 total 6.581421e-07 1.337227e-07 +3 2 2 1 1 total 4.714011e-07 1.140489e-07 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 0.025920 0.002893 -1 1 2 1 1 total 0.022467 0.004039 -2 2 1 1 1 total 0.028922 0.006394 -3 2 2 1 1 total 0.024923 0.009632 +0 1 1 1 1 total 0.020173 0.004383 +1 1 2 1 1 total 0.020397 0.008086 +2 2 1 1 1 total 0.025824 0.003192 +3 2 2 1 1 total 0.020865 0.004879 mesh 1 delayedgroup group in nuclide mean std. dev. x y z -0 1 1 1 1 1 total 0.000004 4.432732e-07 -1 1 1 1 2 1 total 0.000026 2.653319e-06 -2 1 1 1 3 1 total 0.000024 2.402270e-06 -3 1 1 1 4 1 total 0.000054 5.464055e-06 -4 1 1 1 5 1 total 0.000026 2.663025e-06 -5 1 1 1 6 1 total 0.000010 1.038005e-06 -6 1 2 1 1 1 total 0.000004 6.987770e-07 -7 1 2 1 2 1 total 0.000023 4.115234e-06 -8 1 2 1 3 1 total 0.000021 3.816392e-06 -9 1 2 1 4 1 total 0.000049 8.885822e-06 -10 1 2 1 5 1 total 0.000024 4.378290e-06 -11 1 2 1 6 1 total 0.000009 1.745695e-06 -12 2 1 1 1 1 total 0.000005 1.098837e-06 -13 2 1 1 2 1 total 0.000029 6.436855e-06 -14 2 1 1 3 1 total 0.000027 5.926286e-06 -15 2 1 1 4 1 total 0.000061 1.359391e-05 -16 2 1 1 5 1 total 0.000029 6.489015e-06 -17 2 1 1 6 1 total 0.000011 2.574270e-06 -18 2 2 1 1 1 total 0.000004 1.660497e-06 -19 2 2 1 2 1 total 0.000025 9.701974e-06 -20 2 2 1 3 1 total 0.000023 9.005217e-06 -21 2 2 1 4 1 total 0.000054 2.084107e-05 -22 2 2 1 5 1 total 0.000026 9.981045e-06 -23 2 2 1 6 1 total 0.000010 3.987979e-06 +0 1 1 1 1 1 total 0.000005 1.004627e-06 +1 1 1 1 2 1 total 0.000025 5.397916e-06 +2 1 1 1 3 1 total 0.000025 5.274991e-06 +3 1 1 1 4 1 total 0.000058 1.230430e-05 +4 1 1 1 5 1 total 0.000026 5.548686e-06 +5 1 1 1 6 1 total 0.000011 2.306996e-06 +6 1 2 1 1 1 total 0.000004 1.779307e-06 +7 1 2 1 2 1 total 0.000024 9.648451e-06 +8 1 2 1 3 1 total 0.000024 9.479454e-06 +9 1 2 1 4 1 total 0.000056 2.231366e-05 +10 1 2 1 5 1 total 0.000026 1.028267e-05 +11 1 2 1 6 1 total 0.000011 4.268087e-06 +12 2 1 1 1 1 total 0.000006 7.462931e-07 +13 2 1 1 2 1 total 0.000031 3.842628e-06 +14 2 1 1 3 1 total 0.000031 3.666739e-06 +15 2 1 1 4 1 total 0.000071 8.228956e-06 +16 2 1 1 5 1 total 0.000031 3.416811e-06 +17 2 1 1 6 1 total 0.000013 1.429034e-06 +18 2 2 1 1 1 total 0.000005 1.049553e-06 +19 2 2 1 2 1 total 0.000025 5.392757e-06 +20 2 2 1 3 1 total 0.000024 5.145263e-06 +21 2 2 1 4 1 total 0.000056 1.156820e-05 +22 2 2 1 5 1 total 0.000025 4.880083e-06 +23 2 2 1 6 1 total 0.000010 2.037276e-06 mesh 1 delayedgroup group out nuclide mean std. dev. x y z 0 1 1 1 1 1 total 0.0 0.000000 1 1 1 1 2 1 total 0.0 0.000000 2 1 1 1 3 1 total 0.0 0.000000 -3 1 1 1 4 1 total 1.0 1.414214 -4 1 1 1 5 1 total 0.0 0.000000 +3 1 1 1 4 1 total 0.0 0.000000 +4 1 1 1 5 1 total 1.0 1.414214 5 1 1 1 6 1 total 0.0 0.000000 6 1 2 1 1 1 total 0.0 0.000000 7 1 2 1 2 1 total 0.0 0.000000 @@ -184,53 +184,53 @@ 23 2 2 1 6 1 total 0.0 0.000000 mesh 1 delayedgroup group in nuclide mean std. dev. x y z -0 1 1 1 1 1 total 0.000166 0.000023 -1 1 1 1 2 1 total 0.000989 0.000136 -2 1 1 1 3 1 total 0.000907 0.000123 -3 1 1 1 4 1 total 0.002087 0.000282 -4 1 1 1 5 1 total 0.001014 0.000137 -5 1 1 1 6 1 total 0.000400 0.000054 -6 1 2 1 1 1 total 0.000167 0.000030 -7 1 2 1 2 1 total 0.001002 0.000178 -8 1 2 1 3 1 total 0.000926 0.000165 -9 1 2 1 4 1 total 0.002149 0.000384 -10 1 2 1 5 1 total 0.001056 0.000189 -11 1 2 1 6 1 total 0.000417 0.000076 -12 2 1 1 1 1 total 0.000171 0.000039 -13 2 1 1 2 1 total 0.001003 0.000226 -14 2 1 1 3 1 total 0.000918 0.000208 -15 2 1 1 4 1 total 0.002100 0.000477 -16 2 1 1 5 1 total 0.000996 0.000228 -17 2 1 1 6 1 total 0.000394 0.000090 -18 2 2 1 1 1 total 0.000171 0.000082 -19 2 2 1 2 1 total 0.001007 0.000480 -20 2 2 1 3 1 total 0.000929 0.000445 -21 2 2 1 4 1 total 0.002143 0.001028 -22 2 2 1 5 1 total 0.001026 0.000492 -23 2 2 1 6 1 total 0.000408 0.000196 +0 1 1 1 1 1 total 0.000227 0.000061 +1 1 1 1 2 1 total 0.001228 0.000327 +2 1 1 1 3 1 total 0.001206 0.000320 +3 1 1 1 4 1 total 0.002835 0.000748 +4 1 1 1 5 1 total 0.001300 0.000340 +5 1 1 1 6 1 total 0.000540 0.000141 +6 1 2 1 1 1 total 0.000218 0.000074 +7 1 2 1 2 1 total 0.001182 0.000400 +8 1 2 1 3 1 total 0.001160 0.000393 +9 1 2 1 4 1 total 0.002726 0.000926 +10 1 2 1 5 1 total 0.001249 0.000428 +11 1 2 1 6 1 total 0.000519 0.000177 +12 2 1 1 1 1 total 0.000226 0.000033 +13 2 1 1 2 1 total 0.001207 0.000174 +14 2 1 1 3 1 total 0.001175 0.000167 +15 2 1 1 4 1 total 0.002721 0.000378 +16 2 1 1 5 1 total 0.001207 0.000160 +17 2 1 1 6 1 total 0.000502 0.000067 +18 2 2 1 1 1 total 0.000220 0.000068 +19 2 2 1 2 1 total 0.001178 0.000355 +20 2 2 1 3 1 total 0.001150 0.000343 +21 2 2 1 4 1 total 0.002675 0.000783 +22 2 2 1 5 1 total 0.001199 0.000341 +23 2 2 1 6 1 total 0.000499 0.000142 mesh 1 delayedgroup group in nuclide mean std. dev. x y z -0 1 1 1 1 1 total 0.00000 0.000000 -1 1 1 1 2 1 total 0.00000 0.000000 -2 1 1 1 3 1 total 0.00000 0.000000 -3 1 1 1 4 1 total 0.30278 0.428196 -4 1 1 1 5 1 total 0.00000 0.000000 -5 1 1 1 6 1 total 0.00000 0.000000 -6 1 2 1 1 1 total 0.00000 0.000000 -7 1 2 1 2 1 total 0.00000 0.000000 -8 1 2 1 3 1 total 0.00000 0.000000 -9 1 2 1 4 1 total 0.00000 0.000000 -10 1 2 1 5 1 total 0.00000 0.000000 -11 1 2 1 6 1 total 0.00000 0.000000 -12 2 1 1 1 1 total 0.00000 0.000000 -13 2 1 1 2 1 total 0.00000 0.000000 -14 2 1 1 3 1 total 0.00000 0.000000 -15 2 1 1 4 1 total 0.00000 0.000000 -16 2 1 1 5 1 total 0.00000 0.000000 -17 2 1 1 6 1 total 0.00000 0.000000 -18 2 2 1 1 1 total 0.00000 0.000000 -19 2 2 1 2 1 total 0.00000 0.000000 -20 2 2 1 3 1 total 0.00000 0.000000 -21 2 2 1 4 1 total 0.00000 0.000000 -22 2 2 1 5 1 total 0.00000 0.000000 -23 2 2 1 6 1 total 0.00000 0.000000 +0 1 1 1 1 1 total 0.00000 0.00000 +1 1 1 1 2 1 total 0.00000 0.00000 +2 1 1 1 3 1 total 0.00000 0.00000 +3 1 1 1 4 1 total 0.00000 0.00000 +4 1 1 1 5 1 total 0.84949 1.20136 +5 1 1 1 6 1 total 0.00000 0.00000 +6 1 2 1 1 1 total 0.00000 0.00000 +7 1 2 1 2 1 total 0.00000 0.00000 +8 1 2 1 3 1 total 0.00000 0.00000 +9 1 2 1 4 1 total 0.00000 0.00000 +10 1 2 1 5 1 total 0.00000 0.00000 +11 1 2 1 6 1 total 0.00000 0.00000 +12 2 1 1 1 1 total 0.00000 0.00000 +13 2 1 1 2 1 total 0.00000 0.00000 +14 2 1 1 3 1 total 0.00000 0.00000 +15 2 1 1 4 1 total 0.00000 0.00000 +16 2 1 1 5 1 total 0.00000 0.00000 +17 2 1 1 6 1 total 0.00000 0.00000 +18 2 2 1 1 1 total 0.00000 0.00000 +19 2 2 1 2 1 total 0.00000 0.00000 +20 2 2 1 3 1 total 0.00000 0.00000 +21 2 2 1 4 1 total 0.00000 0.00000 +22 2 2 1 5 1 total 0.00000 0.00000 +23 2 2 1 6 1 total 0.00000 0.00000 diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index 5f4bebb0d..d716a171c 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -1 +1 @@ -96d2f92c6017e62d688e3ca245c2ff7904a92816d188d096829440cd8fcfd4f8639dc67c12d54b324081c19e1c7dd2e79c1caeac53d7826399f11766826d5410 \ No newline at end of file +b112ccd96ff89c706c62b4c334021e352de78fc494c4c0513ca7a924bb2e9b18139c5b4fd50c9a92fff07525cc1d9b1464fcc5ca24293e6fd3a8d03d9e12b171 \ No newline at end of file diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat index 108c7ee80..6b3bbf4d5 100644 --- a/tests/test_tallies/results_true.dat +++ b/tests/test_tallies/results_true.dat @@ -1 +1 @@ -a6afd2f11affce2467d77b8477881ab20091f67df4f632226ec2dd5d4cd7fabb9ac3e182563bb467ed249e4b3fe95b319cb688d653757f8ea154759b8a7f50e1 \ No newline at end of file +1e0ce3915c37809cf147e5c0634b696bafd99fa3f41573235d884e54a9a62e8440234a4ee1bf4ce94b363e7afde2ded7c0fd22d0dd897b8aacf1ee017b974449 \ No newline at end of file diff --git a/tests/test_tally_aggregation/inputs_true.dat b/tests/test_tally_aggregation/inputs_true.dat index 6d2990754..d0f617dbf 100644 --- a/tests/test_tally_aggregation/inputs_true.dat +++ b/tests/test_tally_aggregation/inputs_true.dat @@ -1 +1 @@ -4a4e481b9af3612c71bdc93245011555807061bbd9d9be4c5b399f2c38820d38d9ce3ac6255d045415a216737eac65fb0f0b6e331e49b90bf8edc814d0f2f0e8 \ No newline at end of file +e6604f58a6a3115b902364ac9c73fcb8babb9df0c527b0b04780fefda61fb41856cffac548a1f81fba5fdac8ab55dd00f683b724c732ed54506933e9d943574d \ No newline at end of file diff --git a/tests/test_tally_aggregation/results_true.dat b/tests/test_tally_aggregation/results_true.dat index 6c2d7a519..929fc9fcf 100644 --- a/tests/test_tally_aggregation/results_true.dat +++ b/tests/test_tally_aggregation/results_true.dat @@ -1 +1 @@ -840d2648f9ba782926c71baa84e5a2ad31331e156740a3d1e9d86af8f1f0d301ef8c0f69474975d365dbcf8d229a68c62d3e60286d18045e5254373f4e1010bf \ No newline at end of file +2fba1e497435b63ee277a1c8a0ecc566440c56a25cdfef1b790b839fc93da1414bd1af38a63a64fcf3e868d1347ad138f528ae12c1067aea29733670d2c9128f \ No newline at end of file diff --git a/tests/test_tally_aggregation/test_tally_aggregation.py b/tests/test_tally_aggregation/test_tally_aggregation.py index 955d7fecb..f1915b261 100644 --- a/tests/test_tally_aggregation/test_tally_aggregation.py +++ b/tests/test_tally_aggregation/test_tally_aggregation.py @@ -16,9 +16,9 @@ class TallyAggregationTestHarness(PyAPITestHarness): self._input_set.settings.output = {'summary': True} # Initialize the nuclides + u234 = openmc.Nuclide('U234') u235 = openmc.Nuclide('U235') u238 = openmc.Nuclide('U238') - pu239 = openmc.Nuclide('Pu239') # Initialize the filters energy_filter = openmc.EnergyFilter([0.0, 0.253e-6, 1.0e-3, 1.0, 20.0]) @@ -28,7 +28,7 @@ class TallyAggregationTestHarness(PyAPITestHarness): tally = openmc.Tally(name='distribcell tally') tally.filters = [energy_filter, distrib_filter] tally.scores = ['nu-fission', 'total'] - tally.nuclides = [u235, u238, pu239] + tally.nuclides = [u234, u235, u238] tallies_file = openmc.Tallies([tally]) # Export tallies to file @@ -59,7 +59,7 @@ class TallyAggregationTestHarness(PyAPITestHarness): outstr += ', '.join(map(str, tally_sum.std_dev)) # Sum across all nuclides - tally_sum = tally.summation(nuclides=['U235', 'U238', 'Pu239']) + tally_sum = tally.summation(nuclides=['U234', 'U235', 'U238']) outstr += ', '.join(map(str, tally_sum.mean)) outstr += ', '.join(map(str, tally_sum.std_dev)) diff --git a/tests/test_tally_arithmetic/inputs_true.dat b/tests/test_tally_arithmetic/inputs_true.dat index b0da0ed24..8cf6c4ff5 100644 --- a/tests/test_tally_arithmetic/inputs_true.dat +++ b/tests/test_tally_arithmetic/inputs_true.dat @@ -1 +1 @@ -6747131dad4c1efd8c87d857ce9127cb16ec883fdf5fabe309d82732d47851424273e9ad4878a335555d8b25bb0c0a15e68ac30e355ae346e9885c499e9ec979 \ No newline at end of file +530c54691a4f2df131627849dd58c5e6f331ede786c1ffff0b1126c4c0fe64929a15113699103a05f2c6178ec90be4a3d6fc12b80c32ec608ed4a59c965ef2e4 \ No newline at end of file diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py index 6e52313ef..d28db3e95 100644 --- a/tests/test_tally_arithmetic/test_tally_arithmetic.py +++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py @@ -19,9 +19,9 @@ class TallyArithmeticTestHarness(PyAPITestHarness): tallies_file = openmc.Tallies() # Initialize the nuclides + u234 = openmc.Nuclide('U234') u235 = openmc.Nuclide('U235') u238 = openmc.Nuclide('U238') - pu239 = openmc.Nuclide('Pu239') # Initialize Mesh mesh = openmc.Mesh(mesh_id=1) @@ -40,7 +40,7 @@ class TallyArithmeticTestHarness(PyAPITestHarness): tally = openmc.Tally(name='tally 1') tally.filters = [material_filter, energy_filter, distrib_filter] tally.scores = ['nu-fission', 'total'] - tally.nuclides = [u235, pu239] + tally.nuclides = [u234, u235] tallies_file.append(tally) tally = openmc.Tally(name='tally 2') diff --git a/tests/test_tally_slice_merge/inputs_true.dat b/tests/test_tally_slice_merge/inputs_true.dat index 7a747e6b5..ef2595fbf 100644 --- a/tests/test_tally_slice_merge/inputs_true.dat +++ b/tests/test_tally_slice_merge/inputs_true.dat @@ -1 +1 @@ -144dd4059444fad5e2e4fa20681fbdc74c0e5cbf3265104a0b49d87c768798eaeb25e3c6c795bcac2eebdde784c51588006d62c9f25be96dda5c51011a76b7c1 \ No newline at end of file +41b06951808ee7ed1e44defbaadff7305737ae494e8e98c682aefa272c0682cbd8803dd2658d8e7d78e3f5f7fb1e9f0f5e443a9374a68474bd25d0d649ac6924 \ No newline at end of file diff --git a/tests/test_tally_slice_merge/results_true.dat b/tests/test_tally_slice_merge/results_true.dat index 278d8ee10..0f43b35e1 100644 --- a/tests/test_tally_slice_merge/results_true.dat +++ b/tests/test_tally_slice_merge/results_true.dat @@ -1,36 +1,36 @@ cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 0.00e+00 6.25e-07 U235 fission 1.08e-01 7.94e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 0.00e+00 6.25e-07 U235 nu-fission 2.64e-01 1.94e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 0.00e+00 6.25e-07 U238 fission 1.51e-07 1.00e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 0.00e+00 6.25e-07 U238 nu-fission 3.76e-07 2.50e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 6.25e-07 2.00e+01 U235 fission 3.12e-02 2.56e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 6.25e-07 2.00e+01 U235 nu-fission 7.65e-02 6.24e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 6.25e-07 2.00e+01 U238 fission 2.00e-02 1.30e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 6.25e-07 2.00e+01 U238 nu-fission 5.56e-02 3.78e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 0.00e+00 6.25e-07 U235 fission 4.43e-02 7.21e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 0.00e+00 6.25e-07 U235 nu-fission 1.08e-01 1.76e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 0.00e+00 6.25e-07 U238 fission 6.14e-08 9.64e-09 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 0.00e+00 6.25e-07 U238 nu-fission 1.53e-07 2.40e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 6.25e-07 2.00e+01 U235 fission 1.39e-02 1.06e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 6.25e-07 2.00e+01 U235 nu-fission 3.40e-02 2.61e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 6.25e-07 2.00e+01 U238 fission 9.72e-03 1.21e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 27 6.25e-07 2.00e+01 U238 nu-fission 2.71e-02 3.80e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 21 0.00e+00 6.25e-07 U235 fission 1.08e-01 7.94e-03 -1 21 0.00e+00 6.25e-07 U235 nu-fission 2.64e-01 1.94e-02 -2 21 0.00e+00 6.25e-07 U238 fission 1.51e-07 1.00e-08 -3 21 0.00e+00 6.25e-07 U238 nu-fission 3.76e-07 2.50e-08 -4 21 6.25e-07 2.00e+01 U235 fission 3.12e-02 2.56e-03 -5 21 6.25e-07 2.00e+01 U235 nu-fission 7.65e-02 6.24e-03 -6 21 6.25e-07 2.00e+01 U238 fission 2.00e-02 1.30e-03 -7 21 6.25e-07 2.00e+01 U238 nu-fission 5.56e-02 3.78e-03 -8 27 0.00e+00 6.25e-07 U235 fission 4.43e-02 7.21e-03 -9 27 0.00e+00 6.25e-07 U235 nu-fission 1.08e-01 1.76e-02 -10 27 0.00e+00 6.25e-07 U238 fission 6.14e-08 9.64e-09 -11 27 0.00e+00 6.25e-07 U238 nu-fission 1.53e-07 2.40e-08 -12 27 6.25e-07 2.00e+01 U235 fission 1.39e-02 1.06e-03 -13 27 6.25e-07 2.00e+01 U235 nu-fission 3.40e-02 2.61e-03 -14 27 6.25e-07 2.00e+01 U238 fission 9.72e-03 1.21e-03 -15 27 6.25e-07 2.00e+01 U238 nu-fission 2.71e-02 3.80e-03 +0 21 0.00e+00 6.25e-07 U235 fission 2.36e-01 2.26e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 0.00e+00 6.25e-07 U235 nu-fission 5.74e-01 5.51e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 0.00e+00 6.25e-07 U238 fission 3.19e-07 3.06e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 0.00e+00 6.25e-07 U238 nu-fission 7.96e-07 7.63e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 6.25e-07 2.00e+01 U235 fission 2.82e-02 1.82e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 6.25e-07 2.00e+01 U235 nu-fission 6.92e-02 4.42e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 6.25e-07 2.00e+01 U238 fission 1.98e-02 1.74e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 6.25e-07 2.00e+01 U238 nu-fission 5.54e-02 4.65e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 0.00e+00 6.25e-07 U235 fission 1.11e-01 1.16e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 0.00e+00 6.25e-07 U235 nu-fission 2.71e-01 2.83e-02 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 0.00e+00 6.25e-07 U238 fission 1.53e-07 1.67e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 0.00e+00 6.25e-07 U238 nu-fission 3.81e-07 4.15e-08 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 6.25e-07 2.00e+01 U235 fission 2.00e-02 2.72e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 6.25e-07 2.00e+01 U235 nu-fission 4.89e-02 6.62e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 6.25e-07 2.00e+01 U238 fission 1.06e-02 6.18e-04 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 27 6.25e-07 2.00e+01 U238 nu-fission 2.91e-02 1.85e-03 cell energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 21 0.00e+00 6.25e-07 U235 fission 2.36e-01 2.26e-02 +1 21 0.00e+00 6.25e-07 U235 nu-fission 5.74e-01 5.51e-02 +2 21 0.00e+00 6.25e-07 U238 fission 3.19e-07 3.06e-08 +3 21 0.00e+00 6.25e-07 U238 nu-fission 7.96e-07 7.63e-08 +4 21 6.25e-07 2.00e+01 U235 fission 2.82e-02 1.82e-03 +5 21 6.25e-07 2.00e+01 U235 nu-fission 6.92e-02 4.42e-03 +6 21 6.25e-07 2.00e+01 U238 fission 1.98e-02 1.74e-03 +7 21 6.25e-07 2.00e+01 U238 nu-fission 5.54e-02 4.65e-03 +8 27 0.00e+00 6.25e-07 U235 fission 1.11e-01 1.16e-02 +9 27 0.00e+00 6.25e-07 U235 nu-fission 2.71e-01 2.83e-02 +10 27 0.00e+00 6.25e-07 U238 fission 1.53e-07 1.67e-08 +11 27 0.00e+00 6.25e-07 U238 nu-fission 3.81e-07 4.15e-08 +12 27 6.25e-07 2.00e+01 U235 fission 2.00e-02 2.72e-03 +13 27 6.25e-07 2.00e+01 U235 nu-fission 4.89e-02 6.62e-03 +14 27 6.25e-07 2.00e+01 U238 fission 1.06e-02 6.18e-04 +15 27 6.25e-07 2.00e+01 U238 nu-fission 2.91e-02 1.85e-03 sum(distribcell) energy low [MeV] energy high [MeV] nuclide score mean std. dev. 0 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U235 fission 0.00e+00 0.00e+00 1 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U235 nu-fission 0.00e+00 0.00e+00 @@ -49,19 +49,19 @@ 14 (500, 5000, 50000) 6.25e-07 2.00e+01 U238 fission 0.00e+00 0.00e+00 15 (500, 5000, 50000) 6.25e-07 2.00e+01 U238 nu-fission 0.00e+00 0.00e+00 sum(mesh) energy low [MeV] energy high [MeV] nuclide score mean std. dev. -0 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U235 fission 8.54e-03 1.30e-03 -1 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U235 nu-fission 2.08e-02 3.17e-03 -2 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U238 fission 1.21e-08 1.74e-09 -3 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U238 nu-fission 3.01e-08 4.34e-09 -4 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U235 fission 2.20e-03 6.05e-04 -5 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U235 nu-fission 5.38e-03 1.48e-03 -6 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U238 fission 1.40e-03 7.17e-04 -7 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U238 nu-fission 3.84e-03 1.97e-03 -8 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U235 fission 9.40e-03 1.62e-03 -9 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U235 nu-fission 2.29e-02 3.95e-03 -10 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U238 fission 1.34e-08 2.08e-09 -11 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U238 nu-fission 3.33e-08 5.18e-09 -12 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U235 fission 9.41e-04 2.52e-04 -13 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U235 nu-fission 2.31e-03 6.13e-04 -14 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U238 fission 7.54e-04 3.45e-04 -15 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U238 nu-fission 2.12e-03 1.02e-03 +0 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U235 fission 1.12e-02 2.45e-03 +1 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U235 nu-fission 2.74e-02 5.96e-03 +2 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U238 fission 1.56e-08 3.04e-09 +3 ((1, 1, 1), (1, 2, 1)) 0.00e+00 6.25e-07 U238 nu-fission 3.89e-08 7.57e-09 +4 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U235 fission 1.15e-03 4.71e-04 +5 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U235 nu-fission 2.81e-03 1.15e-03 +6 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U238 fission 3.25e-04 9.99e-05 +7 ((1, 1, 1), (1, 2, 1)) 6.25e-07 2.00e+01 U238 nu-fission 8.81e-04 2.68e-04 +8 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U235 fission 1.44e-02 5.57e-03 +9 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U235 nu-fission 3.51e-02 1.36e-02 +10 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U238 fission 1.95e-08 7.37e-09 +11 ((2, 1, 1), (2, 2, 1)) 0.00e+00 6.25e-07 U238 nu-fission 4.86e-08 1.84e-08 +12 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U235 fission 2.55e-03 1.30e-03 +13 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U235 nu-fission 6.23e-03 3.16e-03 +14 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U238 fission 8.62e-04 4.53e-04 +15 ((2, 1, 1), (2, 2, 1)) 6.25e-07 2.00e+01 U238 nu-fission 2.30e-03 1.21e-03 From b3ec0d345ff4dba107b9a27ca8dd4468c8475316 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 10 Oct 2016 18:11:48 -0400 Subject: [PATCH 24/32] Remove bad English --- src/nuclide_header.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 93832869c..3533fc969 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -194,8 +194,6 @@ module nuclide_header integer, intent(inout) :: method real(8), intent(in) :: tolerance logical, intent(in) :: master ! if this is the master proc - ! (can't use global cuz - ! circular dependance) integer :: i integer :: storage_type From 6976960d17467d37e3c804999de5d273b4a43976 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Oct 2016 14:28:58 -0500 Subject: [PATCH 25/32] Disallow export_to_hdf5() when data came from ENDF --- openmc/data/neutron.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 166ffe42a..c125e3e04 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -377,7 +377,7 @@ class IncidentNeutron(EqualityMixin): return mts def export_to_hdf5(self, path, mode='a'): - """Export table to an HDF5 file. + """Export incident neutron data to an HDF5 file. Parameters ---------- @@ -388,6 +388,10 @@ class IncidentNeutron(EqualityMixin): to the :class:`h5py.File` constructor. """ + # If data come from ENDF, don't allow exporting to HDF5 + if hasattr(self, '_evaluation'): + raise NotImplementedError('Cannot export incident neutron data that ' + 'originated from an ENDF file.') f = h5py.File(path, mode, libver='latest') From 37d230e65bc1621953baa76c90222df7b1376af1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Oct 2016 21:09:50 -0500 Subject: [PATCH 26/32] Fix typo in units description --- openmc/data/fission_energy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 687cb5703..83f61a1ae 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -24,7 +24,7 @@ def _extract_458_data(ev, units='eV'): ev : openmc.data.Evaluation ENDF evaluation units : {'eV', 'MeV'} - The units are used in values returned. + The units that are used in values returned. Returns ------- From 965a1acf60e0d0d4addac6029917fdcbde50db51 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 12 Oct 2016 07:45:16 -0500 Subject: [PATCH 27/32] Make sure MF=6 distributions get applied to fission products correctly --- openmc/data/reaction.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index af903fe8a..922e13ed4 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1025,7 +1025,12 @@ class Reaction(EqualityMixin): if (6, mt) in ev.section: # Product angle-energy distribution - rx.products = _get_products(ev, mt) + for product in _get_products(ev, mt): + if mt in (18, 19, 20, 21, 38) and product.particle == 'neutron': + rx.products[0].applicability = product.applicability + rx.products[0].distribution = product.distribution + else: + rx.products.append(product) elif (4, mt) in ev.section or (5, mt) in ev.section: # Uncorrelated angle-energy distribution From 04e670c71563f868975ea07fdb5d14787466b5c4 Mon Sep 17 00:00:00 2001 From: Qingming He <906459647@qq.com> Date: Sat, 8 Oct 2016 16:53:14 -0400 Subject: [PATCH 28/32] add a create_fission_neutrons option Whether create fission neutrons --- src/global.F90 | 3 +++ src/input_xml.F90 | 13 +++++++++++++ src/physics.F90 | 5 +++-- src/physics_mg.F90 | 2 +- 4 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index f00a1d368..a5720464d 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -352,6 +352,9 @@ module global ! Write out initial source logical :: write_initial_source = .false. + ! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE + logical :: create_fission_neutrons = .true. + ! ============================================================================ ! CMFD VARIABLES diff --git a/src/input_xml.F90 b/src/input_xml.F90 index db08eedaa..c92ceee8c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1102,6 +1102,19 @@ contains end select end if + ! Check whether create fission sites + if (run_mode == MODE_FIXEDSOURCE) then + if (check_for_node(doc, "create_fission_neutrons")) then + call get_node_value(doc, "create_fission_neutrons", temp_str) + temp_str = to_lower(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') then + create_fission_neutrons = .true. + else if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') then + create_fission_neutrons = .false. + end if + end if + end if + ! Close settings XML file call close_xmldoc(doc) diff --git a/src/physics.F90 b/src/physics.F90 index 6a6a9787d..232effebf 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -96,10 +96,11 @@ contains ! absorption (including fission) if (nuc % fissionable) then - call sample_fission(i_nuclide, i_reaction) if (run_mode == MODE_EIGENVALUE) then + call sample_fission(i_nuclide, i_reaction) call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank) - elseif (run_mode == MODE_FIXEDSOURCE) then + elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then + call sample_fission(i_nuclide, i_reaction) call create_fission_sites(p, i_nuclide, i_reaction, & p % secondary_bank, p % n_secondary) end if diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 2e5e467c1..c60fc6e2e 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -73,7 +73,7 @@ contains if (mat % fissionable) then if (run_mode == MODE_EIGENVALUE) then call create_fission_sites(p, fission_bank, n_bank) - elseif (run_mode == MODE_FIXEDSOURCE) then + elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then call create_fission_sites(p, p % secondary_bank, p % n_secondary) end if end if From f42a324143fa6b58372552be6fd62e7555b07282 Mon Sep 17 00:00:00 2001 From: Qingming He <906459647@qq.com> Date: Sat, 8 Oct 2016 16:54:32 -0400 Subject: [PATCH 29/32] support create_fission_neutrons option in Python API --- openmc/settings.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/openmc/settings.py b/openmc/settings.py index 335de839a..b453215ea 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -141,6 +141,8 @@ class Settings(object): The elastic scattering model to use for resonant isotopes volume_calculations : VolumeCalculation or iterable of VolumeCalculation Stochastic volume calculation specifications + create_fission_neutrons : bool + Indicate whether fission neutrons should be created or not. """ @@ -227,6 +229,8 @@ class Settings(object): self._volume_calculations = cv.CheckedList( VolumeCalculation, 'volume calculations') + self._create_fission_neutrons = None + @property def run_mode(self): return self._run_mode @@ -427,6 +431,10 @@ class Settings(object): def volume_calculations(self): return self._volume_calculations + @property + def create_fission_neutrons(self): + return self._create_fission_neutrons + @run_mode.setter def run_mode(self, run_mode): if run_mode not in ['eigenvalue', 'fixed source']: @@ -826,6 +834,12 @@ class Settings(object): self._volume_calculations = cv.CheckedList( VolumeCalculation, 'stochastic volume calculations', vol_calcs) + @create_fission_neutrons.setter + def create_fission_neutrons(self, create_fission_neutrons): + cv.check_type('Whether create fission neutrons', + create_fission_neutrons, bool) + self._create_fission_neutrons = create_fission_neutrons + def _create_run_mode_subelement(self): if self.run_mode == 'eigenvalue': @@ -1121,6 +1135,11 @@ class Settings(object): for r in self.resonance_scattering: elem.append(r.to_xml_element()) + def _create_create_fission_neutrons_subelement(self): + if self._create_fission_neutrons is not None: + elem = ET.SubElement(self._settings_file, "create_fission_neutrons") + elem.text = str(self._create_fission_neutrons).lower() + def export_to_xml(self, path='settings.xml'): """Export simulation settings to an XML file. @@ -1164,6 +1183,7 @@ class Settings(object): self._create_dd_subelement() self._create_resonance_scattering_subelement() self._create_volume_calcs_subelement() + self._create_create_fission_neutrons_subelement() # Clean the indentation in the file to be user-readable clean_xml_indentation(self._settings_file) From f4ac440407fd64d5ec02cdc5bb9364dde55e3fb2 Mon Sep 17 00:00:00 2001 From: Qingming He <906459647@qq.com> Date: Sat, 8 Oct 2016 16:55:29 -0400 Subject: [PATCH 30/32] add a create_fission_neutrons unit test --- .../inputs_true.dat | 1 + .../results_true.dat | 3 + .../test_create_fission_neutrons.py | 80 +++++++++++++++++++ 3 files changed, 84 insertions(+) create mode 100644 tests/test_create_fission_neutrons/inputs_true.dat create mode 100644 tests/test_create_fission_neutrons/results_true.dat create mode 100755 tests/test_create_fission_neutrons/test_create_fission_neutrons.py diff --git a/tests/test_create_fission_neutrons/inputs_true.dat b/tests/test_create_fission_neutrons/inputs_true.dat new file mode 100644 index 000000000..8a1a17b9e --- /dev/null +++ b/tests/test_create_fission_neutrons/inputs_true.dat @@ -0,0 +1 @@ +c3581501d1486293c255390251d20584431a198fbb7c07dc2879dc95dc96e26786d3913ce0db8007f542468fd137ff3267824844c03b4687fdad81ea568c92d7 \ No newline at end of file diff --git a/tests/test_create_fission_neutrons/results_true.dat b/tests/test_create_fission_neutrons/results_true.dat new file mode 100644 index 000000000..be7c5c438 --- /dev/null +++ b/tests/test_create_fission_neutrons/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +sum = 2.056839E+02 +sum_sq = 4.244628E+03 diff --git a/tests/test_create_fission_neutrons/test_create_fission_neutrons.py b/tests/test_create_fission_neutrons/test_create_fission_neutrons.py new file mode 100755 index 000000000..80da831cc --- /dev/null +++ b/tests/test_create_fission_neutrons/test_create_fission_neutrons.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + + +class CreateFissionNeutronsTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Material is composed of H-1 and U-235 + mat = openmc.Material(material_id=1, name='mat') + mat.set_density('atom/b-cm', 0.069335) + mat.add_nuclide('H1', 40.0) + mat.add_nuclide('U235', 1.0) + materials_file = openmc.Materials([mat]) + materials_file.export_to_xml() + + # Cell is box with reflective boundary + x1 = openmc.XPlane(surface_id=1, x0=-1) + x2 = openmc.XPlane(surface_id=2, x0=1) + y1 = openmc.YPlane(surface_id=3, y0=-1) + y2 = openmc.YPlane(surface_id=4, y0=1) + z1 = openmc.ZPlane(surface_id=5, z0=-1) + z2 = openmc.ZPlane(surface_id=6, z0=1) + for surface in [x1, x2, y1, y2, z1, z2]: + surface.boundary_type = 'reflective' + box = openmc.Cell(cell_id=1, name='box') + box.region = +x1 & -x2 & +y1 & -y2 & +z1 & -z2 + box.fill = mat + root = openmc.Universe(universe_id=0, name='root universe') + root.add_cell(box) + geometry = openmc.Geometry(root) + geometry.export_to_xml() + + # Set the running parameters + settings_file = openmc.Settings() + settings_file.run_mode = 'fixed source' + settings_file.batches = 10 + settings_file.particles = 100 + settings_file.create_fission_neutrons = False + bounds = [-1, -1, -1, 1, 1, 1] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + watt_dist = openmc.stats.Watt() + settings_file.source = openmc.source.Source(space=uniform_dist, + energy=watt_dist) + settings_file.export_to_xml() + + # Create tallies + tallies = openmc.Tallies() + tally = openmc.Tally(1) + tally.scores = ['flux'] + tallies.append(tally) + tallies.export_to_xml() + + def _get_results(self): + """Digest info in the statepoint and return as a string.""" + # Read the statepoint file. + sp = openmc.StatePoint(self._sp_name) + + # Write out tally data. + outstr = '' + t = sp.get_tally() + outstr += 'tally {0}:\n'.format(t.id) + outstr += 'sum = {0:12.6E}\n'.format(t.sum[0, 0, 0]) + outstr += 'sum_sq = {0:12.6E}\n'.format(t.sum_sq[0, 0, 0]) + + return outstr + + def _cleanup(self): + super(CreateFissionNeutronsTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): + os.remove(f) + + +if __name__ == '__main__': + harness = CreateFissionNeutronsTestHarness('statepoint.10.h5', True) + harness.main() From 4e1046585f7277faf8b5de7c1c2bb3bfa0a10cba Mon Sep 17 00:00:00 2001 From: Qingming He <906459647@qq.com> Date: Mon, 10 Oct 2016 22:46:18 -0400 Subject: [PATCH 31/32] add docs for --- docs/source/usersguide/input.rst | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 847dc66ff..4eb83d703 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -872,6 +872,19 @@ displayed. This element takes the following attributes: *Default*: 5 +```` Element +------------------------------------- + +The ```` element indicates whether fission neutrons +should be created or not. If this element is set to "true", fission neutrons +will be created; otherwise the fission is treated as capture and no fission +neutron will be created. Note that this option is only applied to fixed source +calculation. For eigenvalue calculation, fission will always be treated as real +fission. + + *Default*: true + + ```` Element ------------------------- From 83dcfb40f877695012edbe1703983699833ea387 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 12 Oct 2016 18:15:03 -0400 Subject: [PATCH 32/32] Fix mixed-unit error in fission Q values --- openmc/data/fission_energy.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 83f61a1ae..924fe0086 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -40,6 +40,9 @@ def _extract_458_data(ev, units='eV'): caution. """ + cv.check_type('evaluation', ev, Evaluation) + cv.check_value('energy units', units, ('eV', 'MeV')) + if not ev.target['fissionable']: # This nuclide isn't fissionable. return None @@ -356,7 +359,7 @@ class FissionEnergyRelease(EqualityMixin): self._neutrinos = energy_release @classmethod - def _from_dictionary(cls, energy_release, incident_neutron): + def _from_dictionary(cls, energy_release, incident_neutron, units='eV'): """Generate fission energy release data from a dictionary. Parameters @@ -368,6 +371,8 @@ class FissionEnergyRelease(EqualityMixin): data, more for polynomial data. incident_neutron : openmc.data.IncidentNeutron Corresponding incident neutron dataset + units : {'eV', 'MeV'} + The energy units used in the returned object. Returns ------- @@ -375,6 +380,8 @@ class FissionEnergyRelease(EqualityMixin): Fission energy release data """ + cv.check_value('energy units', units, ('eV', 'MeV')) + out = cls() # How many coefficients are given for each component? If we only find @@ -425,24 +432,28 @@ class FissionEnergyRelease(EqualityMixin): raise ValueError('Ambiguous prompt/total nu value.') nu = nu[0] + if units == 'eV': + nu_const = 8.07e6 + else: + nu_const = 8.07 if isinstance(nu, Tabulated1D): ENP = deepcopy(nu) ENP.y = (energy_release['ENP'] + 1.307 * nu.x - - 8.07 * (nu.y - nu.y[0])) + - nu_const * (nu.y - nu.y[0])) elif isinstance(nu, Polynomial): if len(nu) == 1: ENP = Polynomial([energy_release['ENP'][0], 1.307]) else: ENP = Polynomial( - [energy_release['ENP'][0], 1.307 - 8.07*nu.coef[1]] + - [-8.07*c for c in nu.coef[2:]]) + [energy_release['ENP'][0], 1.307 - nu_const*nu.coef[1]] + + [-nu_const*c for c in nu.coef[2:]]) out.prompt_neutrons = ENP return out @classmethod - def from_endf(cls, ev, incident_neutron): + def from_endf(cls, ev, incident_neutron, units='eV'): """Generate fission energy release data from an ENDF file. Parameters @@ -451,6 +462,8 @@ class FissionEnergyRelease(EqualityMixin): ENDF evaluation incident_neutron : openmc.data.IncidentNeutron Corresponding incident neutron dataset + units : {'eV', 'MeV'} + The energy units used in the returned object. Returns ------- @@ -458,6 +471,7 @@ class FissionEnergyRelease(EqualityMixin): Fission energy release data """ + cv.check_type('evaluation', ev, Evaluation) # Check to make sure this ENDF file matches the expected isomer. if ev.target['atomic_number'] != incident_neutron.atomic_number: @@ -473,10 +487,10 @@ class FissionEnergyRelease(EqualityMixin): raise ValueError('The ENDF evaluation is not fissionable.') # Read the 458 data from the ENDF file. - value, uncertainty = _extract_458_data(ev) + value, uncertainty = _extract_458_data(ev, units) # Build the object. - return cls._from_dictionary(value, incident_neutron) + return cls._from_dictionary(value, incident_neutron, units) @classmethod def from_hdf5(cls, group):