Merge pull request #1322 from drewejohnson/heating-less-fission

Process non-fission heating from njoy
This commit is contained in:
Paul Romano 2019-09-05 16:33:18 -05:00 committed by GitHub
commit f9f5187c04
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
5 changed files with 130 additions and 39 deletions

View file

@ -1,11 +1,9 @@
import sys
from collections import OrderedDict
from collections.abc import Iterable, Mapping, MutableMapping
from collections.abc import Mapping, MutableMapping
from io import StringIO
from math import log10
from numbers import Integral, Real
import os
import shutil
import tempfile
from warnings import warn
@ -15,7 +13,8 @@ import h5py
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Library, Table, get_table, get_metadata
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record
from .endf import (
Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations)
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D, Sum, ResonancesWithBackground
from .grid import linearize, thin
@ -453,7 +452,7 @@ class IncidentNeutron(EqualityMixin):
if rx.redundant:
photon_rx = any(p.particle == 'photon' for p in rx.products)
keep_mts = (4, 16, 103, 104, 105, 106, 107,
203, 204, 205, 206, 207, 301, 444)
203, 204, 205, 206, 207, 301, 318, 444)
if not (photon_rx or rx.mt in keep_mts):
continue
@ -554,6 +553,20 @@ class IncidentNeutron(EqualityMixin):
fer_group = group['fission_energy_release']
data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
# Rebuild non-fission heating
total_heating = data.reactions.get(301)
fission_heating = data.reactions.get(318)
if total_heating is not None and fission_heating is not None:
non_fission_heating = Reaction(999)
non_fission_heating.redundant = True
for strT, total in total_heating.xs.items():
fission = fission_heating.xs.get(strT)
if fission is None:
continue
non_fission_heating.xs[strT] = Tabulated1D(
total.x, total.y - fission(total.x))
data.reactions[999] = non_fission_heating
return data
@classmethod
@ -801,14 +814,14 @@ class IncidentNeutron(EqualityMixin):
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
kwargs.setdefault('pendf', os.path.join(tmpdir, 'pendf'))
kwargs.setdefault("output_dir", tmpdir)
for key in ("acer", "pendf", "heatr", "broadr", "gaspr", "purr"):
kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key))
kwargs['evaluation'] = evaluation
make_ace(filename, temperatures, **kwargs)
# Create instance from ACE tables within library
lib = Library(kwargs['ace'])
lib = Library(kwargs['acer'])
data = cls.from_ace(lib.tables[0])
for table in lib.tables[1:]:
data.add_temperature_from_ace(table)
@ -817,6 +830,29 @@ class IncidentNeutron(EqualityMixin):
ev = evaluation if evaluation is not None else Evaluation(filename)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)
# Add 318 fission heating data from heatr
non_fission_heating = Reaction(999)
non_fission_heating.redundant = True
fission_heating = Reaction(318)
heatr_evals = get_evaluations(kwargs["heatr"])
for heatr in heatr_evals:
temp = "{}K".format(round(heatr.target["temperature"]))
f318 = StringIO(heatr.section[3, 318])
get_head_record(f318)
_params, fission_kerma = get_tab1_record(f318)
fission_heating.xs[temp] = fission_kerma
total_heating_xs = data.reactions[301].xs.get(temp)
if total_heating_xs is None:
continue
non_fission_heating.xs[temp] = Tabulated1D(
fission_kerma.x,
total_heating_xs(fission_kerma.x) - fission_kerma.y,
breakpoints=fission_kerma.breakpoints,
interpolation=fission_kerma.interpolation)
data.reactions[318] = fission_heating
data.reactions[999] = non_fission_heating
# Add 0K elastic scattering cross section
if '0K' not in data.energy:

View file

@ -1,11 +1,10 @@
import argparse
from collections import namedtuple
from io import StringIO
import os
import shutil
from subprocess import Popen, PIPE, STDOUT, CalledProcessError
import sys
import tempfile
from pathlib import Path
from . import endf
@ -215,11 +214,21 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
heatr=False, purr=False, acer=False, stdout=stdout)
def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
error=0.001, broadr=True, heatr=True, gaspr=True, purr=True,
acer=True, evaluation=None, **kwargs):
def make_ace(filename, temperatures=None, acer=True, xsdir=None,
output_dir=None, pendf=False, error=0.001, broadr=True,
heatr=True, gaspr=True, purr=True, evaluation=None,
**kwargs):
"""Generate incident neutron ACE file from an ENDF file
File names can be passed to
``[acer, xsdir, pendf, broadr, heatr, gaspr, purr]``
to specify the exact output for the given module.
Otherwise, the files will be writen to the current directory
or directory specified by ``output_dir``. Default file
names mirror the variable names, e.g. ``heatr`` output
will be written to a file named ``heatr`` unless otherwise
specified.
Parameters
----------
filename : str
@ -227,24 +236,34 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
temperatures : iterable of float, optional
Temperatures in Kelvin to produce ACE files at. If omitted, data is
produced at room temperature (293.6 K).
ace : str, optional
Path of ACE file to write
acer : bool or str, optional
Flag indicating if acer should be run. If a string is give, write the
resulting ``ace`` file to this location. Path of ACE file to write.
Defaults to ``"ace"``
xsdir : str, optional
Path of xsdir file to write
Path of xsdir file to write. Defaults to ``"xsdir"`` in the same
directory as ``acer``
output_dir : str, optional
Directory to write output for requested modules. If not provided
and at least one of ``[pendf, broadr, heatr, gaspr, purr, acer]``
is ``True``, then write output files to current directory. If given,
must be a path to a directory.
pendf : str, optional
Path of pendf file to write. If omitted, the pendf file is not saved.
error : float, optional
Fractional error tolerance for NJOY processing
broadr : bool, optional
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
broadr : bool or str, optional
Indicating whether to Doppler broaden XS when running NJOY. If string,
write the output tape to this file.
heatr : bool or str, optional
Indicating whether to add heating kerma when running NJOY. If string,
write the output tape to this file.
gaspr : bool or str, optional
Indicating whether to add gas production data when running NJOY.
If string, write the output tape to this file.
purr : bool or str, optional
Indicating whether to add probability table when running NJOY.
If string, write the output tape to this file.
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this argument
indicates which evaluation should be used.
@ -255,8 +274,17 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
------
subprocess.CalledProcessError
If the NJOY process returns with a non-zero status
IOError
If ``output_dir`` does not point to a directory
"""
if output_dir is None:
output_dir = Path()
else:
output_dir = Path(output_dir)
if not output_dir.is_dir():
raise IOError("{} is not a directory".format(output_dir))
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
@ -275,8 +303,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
nendf, npendf = 20, 21
tapein = {nendf: filename}
tapeout = {}
if pendf is not None:
tapeout[npendf] = pendf
if pendf:
tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf
# reconr
commands += _TEMPLATE_RECONR
@ -285,6 +313,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
# broadr
if broadr:
nbroadr = nlast + 1
tapeout[nbroadr] = (output_dir / "broadr") if broadr is True else broadr
commands += _TEMPLATE_BROADR
nlast = nbroadr
@ -292,6 +321,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
if heatr:
nheatr_in = nlast
nheatr = nheatr_in + 1
tapeout[nheatr] = (output_dir / "heatr") if heatr is True else heatr
commands += _TEMPLATE_HEATR
nlast = nheatr
@ -299,6 +329,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
if gaspr:
ngaspr_in = nlast
ngaspr = ngaspr_in + 1
tapeout[ngaspr] = (output_dir / "gaspr") if gaspr is True else gaspr
commands += _TEMPLATE_GASPR
nlast = ngaspr
@ -306,6 +337,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
if purr:
npurr_in = nlast
npurr = npurr_in + 1
tapeout[npurr] = (output_dir / "purr") if purr is True else purr
commands += _TEMPLATE_PURR
nlast = npurr
@ -323,19 +355,21 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
commands += _TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = fname.format(ace, temperature)
tapeout[ndir] = fname.format(xsdir, temperature)
tapeout[nace] = fname.format("ace", temperature)
tapeout[ndir] = fname.format("xsdir", temperature)
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
if acer:
with open(ace, 'w') as ace_file, open(xsdir, 'w') as xsdir_file:
ace = (output_dir / "ace") if acer is True else Path(acer)
xsdir = (ace.parent / "xsdir") if xsdir is None else xsdir
with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
for temperature in temperatures:
# Get contents of ACE file
text = open(fname.format(ace, temperature), 'r').read()
text = open(fname.format("ace", temperature), 'r').read()
# If the target is metastable, make sure that ZAID in the ACE file reflects
# this by adding 400
# If the target is metastable, make sure that ZAID in the ACE
# file reflects this by adding 400
if ev.target['isomeric_state'] > 0:
mass_first_digit = int(text[3])
if mass_first_digit <= 2:
@ -345,13 +379,13 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
ace_file.write(text)
# Concatenate into destination xsdir file
text = open(fname.format(xsdir, temperature), 'r').read()
text = open(fname.format("xsdir", temperature), 'r').read()
xsdir_file.write(text)
# Remove ACE/xsdir files for each temperature
for temperature in temperatures:
os.remove(fname.format(ace, temperature))
os.remove(fname.format(xsdir, temperature))
os.remove(fname.format("ace", temperature))
os.remove(fname.format("xsdir", temperature))
def make_ace_thermal(filename, filename_thermal, temperatures=None,

View file

@ -54,6 +54,7 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)',
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',
318: "fission-heating", 999: "non-fission-heating",
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)})

View file

@ -231,6 +231,10 @@ std::string reaction_name(int mt)
return "(n,Xa)";
} else if (mt == 301) {
return "heating";
} else if (mt == 318) {
return "fission-heating";
} else if (mt == 999) {
return "non-fission-heating";
} else if (mt == 444) {
return "damage-energy";
} else if (mt == COHERENT) {

View file

@ -205,6 +205,22 @@ def test_derived_products(am244):
assert total_neutron.yield_(6e6) == pytest.approx(4.2558)
def test_heating(run_in_tmpdir, am244):
assert 318 in am244.reactions
assert 999 in am244.reactions
# compare values in 999 reaction
total_heating = am244.reactions[301].xs["294K"]
fission_heating = am244.reactions[318].xs["294K"]
expected_y = total_heating.y - fission_heating(total_heating.x)
assert am244.reactions[999].xs["294K"].y == pytest.approx(expected_y)
am244.export_to_hdf5("am244.h5")
read_in = openmc.data.IncidentNeutron.from_hdf5("am244.h5")
assert 318 in read_in.reactions
assert 999 in read_in.reactions
assert read_in.reactions[999].xs["294K"].y == pytest.approx(expected_y)
def test_urr(pu239):
for T, ptable in pu239.urr.items():
assert T.endswith('K')
@ -465,7 +481,7 @@ def test_ace_convert(run_in_tmpdir):
filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-001_H_001.endf')
ace_ascii = 'ace_ascii'
ace_binary = 'ace_binary'
openmc.data.njoy.make_ace(filename, ace=ace_ascii)
openmc.data.njoy.make_ace(filename, acer=ace_ascii)
# Convert to binary
openmc.data.ace.ascii_to_binary(ace_ascii, ace_binary)