Merge pull request #1191 from paulromano/heat-data

Provide gas production, heating, damage energy data in HDF5 files
This commit is contained in:
Jingang Liang 2019-03-12 12:06:11 -04:00 committed by GitHub
commit 830bd4f8e0
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
15 changed files with 187 additions and 72 deletions

View file

@ -200,6 +200,11 @@ Miscellaneous
Multigroup Cross Section Generation
-----------------------------------
- William Boyd, Adam Nelson, Paul K. Romano, Samuel Shaner, Benoit Forget, and
Kord Smith, "`Multigroup Cross-Section Generation with the OpenMC Monte Carlo
Particle Transport Code <https://doi.org/10.1080/00295450.2019.1571828>`_,"
*Nucl. Technol.* (2019).
- William Boyd, Benoit Forget, and Kord Smith, "`A single-step framework to
generate spatially self-shielded multi-group cross sections from Monte Carlo
transport simulations <https://doi.org/10.1016/j.anucene.2018.11.017>`_,"

View file

@ -220,6 +220,16 @@ The following tables show all valid scores:
| |multiplicity from (n,2n), (n,3n), and (n,4n) |
| |reactions. |
+----------------------+---------------------------------------------------+
|H1-production |Total production of H1. |
+----------------------+---------------------------------------------------+
|H2-production |Total production of H2 (deuterium). |
+----------------------+---------------------------------------------------+
|H3-production |Total production of H3 (tritium). |
+----------------------+---------------------------------------------------+
|He3-production |Total production of He3. |
+----------------------+---------------------------------------------------+
|He4-production |Total production of He4 (alpha particles). |
+----------------------+---------------------------------------------------+
.. table:: **Miscellaneous scores: units are indicated for each.**
@ -246,6 +256,10 @@ The following tables show all valid scores:
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
+----------------------+---------------------------------------------------+
|heating |Total neutron heating in units of eV per source |
| |particle. This corresponds to MT=301 produced by |
| |NJOY's HEATR module. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
| |fission product kinetic energy, prompt and delayed |
@ -281,3 +295,7 @@ The following tables show all valid scores:
|decay-rate |The delayed-nu-fission-weighted decay rate where |
| |the decay rate is in units of inverse seconds. |
+----------------------+---------------------------------------------------+
|damage-energy |Damage energy production in units of eV per source |
| |particle. This corresponds to MT=444 produced by |
| |NJOY's HEATR module. |
+----------------------+---------------------------------------------------+

View file

@ -225,6 +225,13 @@ constexpr int N_3P {197};
constexpr int N_N3P {198};
constexpr int N_3N2PA {199};
constexpr int N_5N2P {200};
constexpr int N_XP {203};
constexpr int N_XD {204};
constexpr int N_XT {205};
constexpr int N_X3HE {206};
constexpr int N_XA {207};
constexpr int HEATING {301};
constexpr int DAMAGE_ENERGY {444};
constexpr int COHERENT {502};
constexpr int INCOHERENT {504};
constexpr int PAIR_PROD_ELEC {515};

View file

@ -66,7 +66,7 @@ SUM_RULES = {1: [2, 3],
106: list(range(750, 800)),
107: list(range(800, 850))}
ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]) ?(\d+)')
def float_endf(s):
@ -89,12 +89,15 @@ def float_endf(s):
The number
"""
return float(ENDF_FLOAT_RE.sub(r'\1e\2', s))
return float(ENDF_FLOAT_RE.sub(r'\1e\2\3', s))
def _int_endf(s):
"""Convert string to int. Used for INTG records where blank entries
indicate a 0.
def int_endf(s):
"""Convert string of integer number in ENDF to int.
The ENDF-6 format technically allows integers to be represented by a field
of all blanks. This function acts like int(s) except when s is a string of
all whitespace, in which case zero is returned.
Parameters
----------
@ -106,8 +109,7 @@ def _int_endf(s):
integer
The number or 0
"""
s = s.strip()
return int(s) if s else 0
return 0 if s.isspace() else int(s)
def get_text_record(file_obj):
@ -127,35 +129,35 @@ def get_text_record(file_obj):
return file_obj.readline()[:66]
def get_cont_record(file_obj, skipC=False):
def get_cont_record(file_obj, skip_c=False):
"""Return data from a CONT record in an ENDF-6 file.
Parameters
----------
file_obj : file-like object
ENDF-6 file to read from
skipC : bool
skip_c : bool
Determine whether to skip the first two quantities (C1, C2) of the CONT
record.
Returns
-------
list
tuple
The six items within the CONT record
"""
line = file_obj.readline()
if skipC:
if skip_c:
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]
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
N1 = int_endf(line[44:55])
N2 = int_endf(line[55:66])
return (C1, C2, L1, L2, N1, N2)
def get_head_record(file_obj):
@ -168,18 +170,18 @@ def get_head_record(file_obj):
Returns
-------
list
tuple
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]
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
N1 = int_endf(line[44:55])
N2 = int_endf(line[55:66])
return (ZA, AWR, L1, L2, N1, N2)
def get_list_record(file_obj):
@ -233,10 +235,10 @@ def get_tab1_record(file_obj):
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])
L1 = int_endf(line[22:33])
L2 = int_endf(line[33:44])
n_regions = int_endf(line[44:55])
n_pairs = int_endf(line[55:66])
params = [C1, C2, L1, L2]
# Read the interpolation region data, namely NBT and INT
@ -247,8 +249,8 @@ def get_tab1_record(file_obj):
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])
breakpoints[m] = int_endf(line[0:11])
interpolation[m] = int_endf(line[11:22])
line = line[22:]
m += 1
@ -306,9 +308,9 @@ def get_intg_record(file_obj):
"""
# determine how many items are in list and NDIGIT
items = get_cont_record(file_obj)
ndigit = int(items[2])
npar = int(items[3]) # Number of parameters
nlines = int(items[4]) # Lines to read
ndigit = items[2]
npar = items[3] # Number of parameters
nlines = items[4] # Lines to read
NROW_RULES = {2: 18, 3: 12, 4: 11, 5: 9, 6: 8}
nrow = NROW_RULES[ndigit]
@ -316,13 +318,13 @@ def get_intg_record(file_obj):
corr = np.identity(npar)
for i in range(nlines):
line = file_obj.readline()
ii = _int_endf(line[:5]) - 1 # -1 to account for 0 indexing
jj = _int_endf(line[5:10]) - 1
ii = int_endf(line[:5]) - 1 # -1 to account for 0 indexing
jj = int_endf(line[5:10]) - 1
factor = 10**ndigit
for j in range(nrow):
if jj+j >= ii:
break
element = _int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)])
element = int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)])
if element > 0:
corr[ii, jj] = (element+0.5)/factor
elif element < 0:
@ -507,16 +509,7 @@ class Evaluation(object):
# File numbers, reaction designations, and number of records
for i in range(NXC):
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
_, _, mf, mt, nc, mod = get_cont_record(file_obj, skip_c=True)
self.reaction_list.append((mf, mt, nc, mod))
@property

View file

@ -448,11 +448,13 @@ class IncidentNeutron(EqualityMixin):
for rx in self.reactions.values():
# Skip writing redundant reaction if it doesn't have photon
# production or is a summed transmutation reaction. MT=4 is also
# sometimes needed for probability tables.
# sometimes needed for probability tables. Also write gas
# production, heating, and damage energy production.
if rx.redundant:
photon_rx = any(p.particle == 'photon' for p in rx.products)
transmutation_rx = (rx.mt in (16, 103, 104, 105, 106, 107))
if not (photon_rx or transmutation_rx or rx.mt == 4):
keep_mts = (4, 16, 103, 104, 105, 106, 107,
203, 204, 205, 206, 207, 301, 444)
if not (photon_rx or rx.mt in keep_mts):
continue
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
@ -615,13 +617,14 @@ class IncidentNeutron(EqualityMixin):
# Read energy grid
n_energy = ace.nxs[3]
energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy]*EV_PER_MEV
i = ace.jxs[1]
energy = ace.xss[i : i + n_energy]*EV_PER_MEV
data.energy[strT] = energy
total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2 * n_energy]
absorption_xs = ace.xss[ace.jxs[1] + 2 * n_energy:ace.jxs[1] +
3 * n_energy]
total_xs = ace.xss[i + n_energy : i + 2*n_energy]
absorption_xs = ace.xss[i + 2*n_energy : i + 3*n_energy]
heating_number = ace.xss[i + 3*n_energy : i + 4*n_energy]*EV_PER_MEV
# Create redundant reactions (total and absorption)
# Create redundant reactions (total, absorption, and heating)
total = Reaction(1)
total.xs[strT] = Tabulated1D(energy, total_xs)
total.redundant = True
@ -633,13 +636,15 @@ class IncidentNeutron(EqualityMixin):
absorption.redundant = True
data.reactions[101] = absorption
heating = Reaction(301)
heating.xs[strT] = Tabulated1D(energy, heating_number)
heating.redundant = True
data.reactions[301] = heating
# Read each reaction
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
rx = Reaction.from_ace(ace, i)
# Don't include gas production / damage cross sections
if 200 < rx.mt < 219 or rx.mt == 444:
continue
data.reactions[rx.mt] = rx
# Some photon production reactions may be assigned to MTs that don't
@ -690,6 +695,8 @@ class IncidentNeutron(EqualityMixin):
mts = data.get_reaction_components(rx.mt)
if mts != [rx.mt]:
rx.redundant = True
if rx.mt in (203, 204, 205, 206, 207, 444):
rx.redundant = True
# Read unresolved resonance probability tables
urr = ProbabilityTables.from_ace(ace)
@ -784,16 +791,19 @@ class IncidentNeutron(EqualityMixin):
return data
@classmethod
def from_njoy(cls, filename, temperatures=None, **kwargs):
def from_njoy(cls, filename, temperatures=None, evaluation=None, **kwargs):
"""Generate incident neutron data by running NJOY.
Parameters
----------
filename : str
Path to ENDF evaluation
Path to ENDF file
temperatures : iterable of float
Temperatures in Kelvin to produce data at. If omitted, data is
produced at room temperature (293.6 K)
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this
argument indicates which evaluation to use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace`
@ -808,6 +818,7 @@ class IncidentNeutron(EqualityMixin):
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
pendf_file = os.path.join(tmpdir, 'pendf')
kwargs['evaluation'] = evaluation
make_ace(filename, temperatures, ace_file, xsdir_file,
pendf_file, **kwargs)
@ -818,7 +829,7 @@ class IncidentNeutron(EqualityMixin):
data.add_temperature_from_ace(table)
# Add fission energy release data
ev = Evaluation(filename)
ev = evaluation if evaluation is not None else Evaluation(filename)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)

View file

@ -72,8 +72,13 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
{mat} 3 /
302 318 402 /
{mat} 4 /
302 318 402 444 /
"""
_TEMPLATE_GASPR = """
gaspr / %%%%%%%%%%%%%%%%%%%%%%%%% Add gas production %%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {ngaspr_in} {ngaspr} /
"""
_TEMPLATE_PURR = """
@ -186,7 +191,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
"""Generate ACE file from an ENDF file
"""Generate pointwise ENDF file from an ENDF file
Parameters
----------
@ -211,8 +216,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
error=0.001, broadr=True, heatr=True, purr=True, acer=True,
**kwargs):
error=0.001, broadr=True, heatr=True, gaspr=True, purr=True,
acer=True, evaluation=None, **kwargs):
"""Generate incident neutron ACE file from an ENDF file
Parameters
@ -234,10 +239,15 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
Indicating whether to Doppler broaden XS when running NJOY
heatr : bool, optional
Indicating whether to add heating kerma when running NJOY
gaspr : bool, optional
Indicating whether to add gas production data when running NJOY
purr : bool, optional
Indicating whether to add probability table when running NJOY
acer : bool, optional
Indicating whether to generate ACE file when running NJOY
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this argument
indicates which evaluation should be used.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -247,7 +257,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
If the NJOY process returns with a non-zero status
"""
ev = endf.Evaluation(filename)
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
@ -285,6 +295,13 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
commands += _TEMPLATE_HEATR
nlast = nheatr
# gaspr
if gaspr:
ngaspr_in = nlast
ngaspr = ngaspr_in + 1
commands += _TEMPLATE_GASPR
nlast = ngaspr
# purr
if purr:
npurr_in = nlast
@ -338,7 +355,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', error=0.001, **kwargs):
ace='ace', xsdir='xsdir', error=0.001, evaluation=None,
evaluation_thermal=None, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
@ -356,6 +374,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
Path of xsdir file to write
error : float, optional
Fractional error tolerance for NJOY processing
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.
evaluation_thermal : openmc.data.endf.Evaluation, optional
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
@ -365,11 +389,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
If the NJOY process returns with a non-zero status
"""
ev = endf.Evaluation(filename)
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
ev_thermal = endf.Evaluation(filename_thermal)
ev_thermal = (evaluation_thermal if evaluation_thermal is not None
else endf.Evaluation(filename_thermal))
mat_thermal = ev_thermal.material
zsymam_thermal = ev_thermal.target['zsymam']

View file

@ -51,7 +51,9 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
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)',
198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 203: '(n,Xp)',
204: '(n,Xd)', 205: '(n,Xt)', 206: '(n,X3He)', 207: '(n,Xa)',
301: 'heating', 444: 'damage-energy',
649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)',
849: '(n,ac)', 891: '(n,2nc)'}
REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)})
@ -988,6 +990,10 @@ class Reaction(EqualityMixin):
# Read reaction cross section
xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy]
# For damage energy production, convert to eV
if mt == 444:
xs *= EV_PER_MEV
# Fix negatives -- known issue for Y89 in JEFF 3.2
if np.any(xs < 0.0):
warn("Negative cross sections found for MT={} in {}. Setting "

View file

@ -610,7 +610,8 @@ class ThermalScattering(EqualityMixin):
return table
@classmethod
def from_njoy(cls, filename, filename_thermal, temperatures=None, **kwargs):
def from_njoy(cls, filename, filename_thermal, temperatures=None,
evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate incident neutron data by running NJOY.
Parameters
@ -623,6 +624,13 @@ class ThermalScattering(EqualityMixin):
Temperatures in Kelvin to produce data at. If omitted, data is
produced at all temperatures in the ENDF thermal scattering
sublibrary.
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.
evaluation_thermal : openmc.data.endf.Evaluation, optional
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to
use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal`
@ -636,6 +644,8 @@ class ThermalScattering(EqualityMixin):
# Run NJOY to create an ACE library
ace_file = os.path.join(tmpdir, 'ace')
xsdir_file = os.path.join(tmpdir, 'xsdir')
kwargs['evaluation'] = evaluation
kwargs['evaluation_thermal'] = evaluation_thermal
make_ace_thermal(filename, filename_thermal, temperatures,
ace_file, xsdir_file, **kwargs)

View file

@ -229,8 +229,10 @@ std::string reaction_name(int mt)
return "(n,X3He)";
} else if (mt == 207) {
return "(n,Xa)";
} else if (mt == 301) {
return "heating";
} else if (mt == 444) {
return "(damage)";
return "damage-energy";
} else if (mt == COHERENT) {
return "coherent scatter";
} else if (mt == INCOHERENT) {

View file

@ -197,6 +197,20 @@ score_str_to_int(std::string score_str)
return N_PT;
if (score_str == "(n,da)")
return N_DA;
if (score_str == "(n,Xp)" || score_str == "H1-production")
return N_XP;
if (score_str == "(n,Xd)" || score_str == "H2-production")
return N_XD;
if (score_str == "(n,Xt)" || score_str == "H3-production")
return N_XT;
if (score_str == "(n,X3He)" || score_str == "He3-production")
return N_X3HE;
if (score_str == "(n,Xa)" || score_str == "He4-production")
return N_XA;
if (score_str == "heating")
return HEATING;
if (score_str == "damage-energy")
return DAMAGE_ENERGY;
// So far we have not identified this score string. Check to see if it is a
// deprecated score.

View file

@ -511,4 +511,7 @@
<scores>total</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="33">
<scores>H1-production H2-production H3-production He3-production He4-production heating damage-energy</scores>
</tally>
</tallies>

View file

@ -1 +1 @@
b5edf87cb58db29aa1c38203141df2f055b88aca675aee9673879b579391412abd232cb3d404f60840ac86b936384050576cccb1e3a1fc5c8290b6e890dda74c
de773a799f84348241ebd65bf7beae8114861d8674b652bfd443d578c146a7d37ca38f2efc5d05124fdbb1cf12829907768351e6b2d6b5f4bd2c121417757507

View file

@ -165,6 +165,10 @@ def test_tallies():
all_nuclide_tallies[3].filters = [mesh_filter]
all_nuclide_tallies[3].nuclides = ['U235']
fusion_tally = Tally()
fusion_tally.scores = ['H1-production', 'H2-production', 'H3-production',
'He3-production', 'He4-production', 'heating', 'damage-energy']
model.tallies += [
azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
cellborn_tally, dg_tally, energy_tally, energyout_tally,
@ -174,5 +178,6 @@ def test_tallies():
model.tallies += score_tallies
model.tallies += flux_tallies
model.tallies += all_nuclide_tallies
model.tallies.append(fusion_tally)
harness.main()

View file

@ -0,0 +1,16 @@
from openmc.data import endf
from pytest import approx
def test_float_endf():
assert endf.float_endf('+3.2146') == approx(3.2146)
assert endf.float_endf('.12345') == approx(0.12345)
assert endf.float_endf('6.022+23') == approx(6.022e23)
assert endf.float_endf('6.022-23') == approx(6.022e-23)
assert endf.float_endf(' +1.01+ 2') == approx(101.0)
assert endf.float_endf(' -1.01- 2') == approx(-0.0101)
def test_int_endf():
assert endf.int_endf(' ') == 0
assert endf.int_endf('+4032') == 4032

View file

@ -3,7 +3,7 @@ set -ex
# Download HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget -q -O - https://anl.box.com/shared/static/9jmb8v2cx6kx03s6mbr2ai0mnvx5j79p.xz | tar -C $HOME -xJ
wget -q -O - https://anl.box.com/shared/static/pzutl4i2717yypv12l78l7fn5nmg6grs.xz | tar -C $HOME -xJ
fi
# Download ENDF/B-VII.1 distribution