from collections import namedtuple
from io import StringIO
import os
import shutil
from subprocess import Popen, PIPE, STDOUT, CalledProcessError
import tempfile
from pathlib import Path
from . import endf
# For a given MAT number, give a name for the ACE table and a list of ZAID
# identifiers. This is based on Appendix C in the ENDF manual.
ThermalTuple = namedtuple('ThermalTuple', ['name', 'zaids', 'nmix'])
_THERMAL_DATA = {
1: ThermalTuple('hh2o', [1001], 1),
2: ThermalTuple('parah', [1001], 1),
3: ThermalTuple('orthoh', [1001], 1),
5: ThermalTuple('hyh2', [1001], 1),
7: ThermalTuple('hzrh', [1001], 1),
8: ThermalTuple('hcah2', [1001], 1),
10: ThermalTuple('hice', [1001], 1),
11: ThermalTuple('dd2o', [1002], 1),
12: ThermalTuple('parad', [1002], 1),
13: ThermalTuple('orthod', [1002], 1),
14: ThermalTuple('dice', [1002], 1),
26: ThermalTuple('be', [4009], 1),
27: ThermalTuple('bebeo', [4009], 1),
28: ThermalTuple('bebe2c', [4009], 1),
30: ThermalTuple('graph', [6000, 6012, 6013], 1),
31: ThermalTuple('grph10', [6000, 6012, 6013], 1),
32: ThermalTuple('grph30', [6000, 6012, 6013], 1),
33: ThermalTuple('lch4', [1001], 1),
34: ThermalTuple('sch4', [1001], 1),
35: ThermalTuple('sch4p2', [1001], 1),
37: ThermalTuple('hch2', [1001], 1),
38: ThermalTuple('mesi00', [1001], 1),
39: ThermalTuple('lucite', [1001], 1),
40: ThermalTuple('benz', [1001, 6000, 6012], 2),
42: ThermalTuple('tol00', [1001], 1),
43: ThermalTuple('sisic', [14028, 14029, 14030], 1),
44: ThermalTuple('csic', [6000, 6012, 6013], 1),
45: ThermalTuple('ouo2', [8016, 8017, 8018], 1),
46: ThermalTuple('obeo', [8016, 8017, 8018], 1),
47: ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3),
48: ThermalTuple('osap00', [92238], 1),
49: ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3),
50: ThermalTuple('oice', [8016, 8017, 8018], 1),
51: ThermalTuple('od2o', [8016, 8017, 8018], 1),
52: ThermalTuple('mg24', [12024], 1),
53: ThermalTuple('al27', [13027], 1),
55: ThermalTuple('yyh2', [39089], 1),
56: ThermalTuple('fe56', [26056], 1),
58: ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1),
59: ThermalTuple('si00', [14028], 1),
60: ThermalTuple('asap00', [13027], 1),
71: ThermalTuple('n-un', [7014, 7015], 1),
72: ThermalTuple('u-un', [92238], 1),
75: ThermalTuple('uuo2', [8016, 8017, 8018], 1),
}
def _get_thermal_data(ev, mat):
"""Return appropriate ThermalTuple, accounting for bugs."""
# JEFF assigns MAT=59 to Ca in CaH2 (which is supposed to be silicon).
if ev.info['library'][0] == 'JEFF':
if ev.material == 59:
if 'CaH2' in ''.join(ev.info['description']):
zaids = [20040, 20042, 20043, 20044, 20046, 20048]
return ThermalTuple('cacah2', zaids, 1)
# Before ENDF/B-VIII.0, crystalline graphite was MAT=31
if ev.info['library'] != ('ENDF/B', 8, 0):
if ev.material == 31:
return _THERMAL_DATA[30]
# ENDF/B incorrectly assigns MAT numbers for UO2
#
# Material | ENDF Manual | VII.0 | VII.1 | VIII.0
# ---------|-------------|-------|-------|-------
# O in UO2 | 45 | 75 | 75 | 75
# U in UO2 | 75 | 76 | 48 | 48
if ev.info['library'][0] == 'ENDF/B':
if ev.material == 75:
return _THERMAL_DATA[45]
version = ev.info['library'][1:]
if version in ((7, 1), (8, 0)) and ev.material == 48:
return _THERMAL_DATA[75]
if version == (7, 0) and ev.material == 76:
return _THERMAL_DATA[75]
# If not a problematic material, use the dictionary as is
return _THERMAL_DATA[mat]
_TEMPLATE_RECONR = """
reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%%
{nendf} {npendf}
'{library} PENDF for {zsymam}'/
{mat} 2/
{error}/ err
'{library}: {zsymam}'/
'Processed by NJOY'/
0/
"""
_TEMPLATE_BROADR = """
broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {npendf} {nbroadr}
{mat} {num_temp} 0 0 0. /
{error}/ errthn
{temps}
0/
"""
_TEMPLATE_HEATR = """
heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr} /
{mat} 4 0 0 0 /
302 318 402 444 /
"""
_TEMPLATE_HEATR_LOCAL = """
heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%%
{nendf} {nheatr_in} {nheatr_local} /
{mat} 4 0 0 1 /
302 318 402 444 /
"""
_TEMPLATE_GASPR = """
gaspr / %%%%%%%%%%%%%%%%%%%%%%%%% Add gas production %%%%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {ngaspr_in} {ngaspr} /
"""
_TEMPLATE_PURR = """
purr / %%%%%%%%%%%%%%%%%%%%%%%% Add probability tables %%%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {npurr_in} {npurr} /
{mat} {num_temp} 1 20 64 /
{temps}
1.e10
0/
"""
_TEMPLATE_ACER = """
acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nacer_in} 0 {nace} {ndir}
1 0 1 .{ext} /
'{library}: {zsymam} at {temperature}'/
{mat} {temperature}
1 1/
/
"""
_THERMAL_TEMPLATE_THERMR = """
thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (free gas) %%%%%%%%%%%%%%%
0 {nthermr1_in} {nthermr1}
0 {mat} 12 {num_temp} 1 0 {iform} 1 221 1/
{temps}
{error} {energy_max}
thermr / %%%%%%%%%%%%%%%% Add thermal scattering data (bound) %%%%%%%%%%%%%%%%%%
{nthermal_endf} {nthermr2_in} {nthermr2}
{mat_thermal} {mat} 16 {num_temp} {inelastic} {elastic} {iform} {natom} 222 1/
{temps}
{error} {energy_max}
"""
_THERMAL_TEMPLATE_ACER = """
acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nthermal_acer_in} 0 {nace} {ndir}
2 0 1 .{ext}/
'{library}: {zsymam_thermal} processed by NJOY'/
{mat} {temperature} '{data.name}' /
{zaids} /
222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/
"""
[docs]def run(commands, tapein, tapeout, input_filename=None, stdout=False,
njoy_exec='njoy'):
"""Run NJOY with given commands
Parameters
----------
commands : str
Input commands for NJOY
tapein : dict
Dictionary mapping tape numbers to paths for any input files
tapeout : dict
Dictionary mapping tape numbers to paths for any output files
input_filename : str, optional
File name to write out NJOY input commands
stdout : bool, optional
Whether to display output when running NJOY
njoy_exec : str, optional
Path to NJOY executable
Raises
------
subprocess.CalledProcessError
If the NJOY process returns with a non-zero status
"""
if input_filename is not None:
with open(str(input_filename), 'w') as f:
f.write(commands)
with tempfile.TemporaryDirectory() as tmpdir:
# Copy evaluations to appropriates 'tapes'
for tape_num, filename in tapein.items():
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
shutil.copy(str(filename), tmpfilename)
# Start up NJOY process
njoy = Popen([njoy_exec], cwd=tmpdir, stdin=PIPE, stdout=PIPE,
stderr=STDOUT, universal_newlines=True)
njoy.stdin.write(commands)
njoy.stdin.flush()
lines = []
while True:
# If process is finished, break loop
line = njoy.stdout.readline()
if not line and njoy.poll() is not None:
break
lines.append(line)
if stdout:
# If user requested output, print to screen
print(line, end='')
# Check for error
if njoy.returncode != 0:
raise CalledProcessError(njoy.returncode, njoy_exec,
''.join(lines))
# Copy output files back to original directory
for tape_num, filename in tapeout.items():
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
if os.path.isfile(tmpfilename):
shutil.move(tmpfilename, str(filename))
[docs]def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
"""Generate pointwise ENDF file from an ENDF file
Parameters
----------
filename : str
Path to ENDF file
pendf : str, optional
Path of pointwise ENDF file to write
error : float, optional
Fractional error tolerance for NJOY processing
stdout : bool
Whether to display NJOY standard output
Raises
------
subprocess.CalledProcessError
If the NJOY process returns with a non-zero status
"""
make_ace(filename, pendf=pendf, error=error, broadr=False,
heatr=False, purr=False, acer=False, stdout=stdout)
[docs]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
Path to ENDF file
temperatures : iterable of float, optional
Temperatures in Kelvin to produce ACE files at. If omitted, data is
produced at room temperature (293.6 K).
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. 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 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.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
Raises
------
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']
# Determine name of library
library = '{}-{}.{}'.format(*ev.info['library'])
if temperatures is None:
temperatures = [293.6]
num_temp = len(temperatures)
temps = ' '.join(str(i) for i in temperatures)
# Create njoy commands by modules
commands = ""
nendf, npendf = 20, 21
tapein = {nendf: filename}
tapeout = {}
if pendf:
tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf
# reconr
commands += _TEMPLATE_RECONR
nlast = npendf
# broadr
if broadr:
nbroadr = nlast + 1
tapeout[nbroadr] = (output_dir / "broadr") if broadr is True else broadr
commands += _TEMPLATE_BROADR
nlast = nbroadr
# heatr
if heatr:
nheatr_in = nlast
nheatr_local = nheatr_in + 1
tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \
else heatr + '_local'
commands += _TEMPLATE_HEATR_LOCAL
nheatr = nheatr_local + 1
tapeout[nheatr] = (output_dir / "heatr") if heatr is True else heatr
commands += _TEMPLATE_HEATR
nlast = nheatr
# gaspr
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
# purr
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
commands = commands.format(**locals())
# acer
if acer:
nacer_in = nlast
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = nacer_in + 1 + 2*i
ndir = nace + 1
ext = '{:02}'.format(i + 1)
commands += _TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
if acer:
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 = (output_dir / "ace_{:.1f}".format(temperature)).read_text()
# 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:
text = text[:3] + str(mass_first_digit + 4) + text[4:]
# Concatenate into destination ACE file
ace_file.write(text)
# Concatenate into destination xsdir file
xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
xsdir_file.write(xsdir_in.read_text())
# Remove ACE/xsdir files for each temperature
for temperature in temperatures:
(output_dir / "ace_{:.1f}".format(temperature)).unlink()
(output_dir / "xsdir_{:.1f}".format(temperature)).unlink()
[docs]def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir=None, output_dir=None, error=0.001,
iwt=2, evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
----------
filename : str
Path to ENDF neutron sublibrary file
filename_thermal : str
Path to ENDF thermal scattering sublibrary file
temperatures : iterable of float, optional
Temperatures in Kelvin to produce data at. If omitted, data is produced
at all temperatures given in the ENDF thermal scattering sublibrary.
ace : str, optional
Path of ACE file to write
xsdir : str, optional
Path of xsdir file to write. Defaults to ``"xsdir"`` in the same
directory as ``ace``
output_dir : str, optional
Directory to write ace and xsdir files. If not provided, then write
output files to current directory. If given, must be a path to a
directory.
error : float, optional
Fractional error tolerance for NJOY processing
iwt : int
`iwt` parameter used in NJOR/ACER card 9
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`
Raises
------
subprocess.CalledProcessError
If the NJOY process returns with a non-zero status
"""
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']
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']
# Determine name, isotopes based on MAT number
data = _get_thermal_data(ev_thermal, mat_thermal)
zaids = ' '.join(str(zaid) for zaid in data.zaids[:3])
# Determine name of library
library = '{}-{}.{}'.format(*ev_thermal.info['library'])
# Determine if thermal elastic is present
if (7, 2) in ev_thermal.section:
elastic = 1
mt_elastic = 223
# Determine whether elastic is incoherent (0) or coherent (1)
file_obj = StringIO(ev_thermal.section[7, 2])
elastic_type = endf.get_head_record(file_obj)[2] - 1
else:
elastic = 0
mt_elastic = 0
elastic_type = 0
# Determine number of principal atoms
file_obj = StringIO(ev_thermal.section[7, 4])
items = endf.get_head_record(file_obj)
items, values = endf.get_list_record(file_obj)
energy_max = values[3]
natom = int(values[5])
# Note that the 'iform' parameter is omitted in NJOY 99. We assume that the
# user is using NJOY 2012 or later.
iform = 0
inelastic = 2
# Determine temperatures from MF=7, MT=4 if none were specified
if temperatures is None:
file_obj = StringIO(ev_thermal.section[7, 4])
endf.get_head_record(file_obj)
endf.get_list_record(file_obj)
endf.get_tab2_record(file_obj)
params = endf.get_tab1_record(file_obj)[0]
temperatures = [params[0]]
for i in range(params[2]):
temperatures.append(endf.get_list_record(file_obj)[0][0])
num_temp = len(temperatures)
temps = ' '.join(str(i) for i in temperatures)
# Create njoy commands by modules
commands = ""
nendf, nthermal_endf, npendf = 20, 21, 22
tapein = {nendf: filename, nthermal_endf: filename_thermal}
tapeout = {}
# reconr
commands += _TEMPLATE_RECONR
nlast = npendf
# broadr
nbroadr = nlast + 1
commands += _TEMPLATE_BROADR
nlast = nbroadr
# thermr
nthermr1_in = nlast
nthermr1 = nthermr1_in + 1
nthermr2_in = nthermr1
nthermr2 = nthermr2_in + 1
commands += _THERMAL_TEMPLATE_THERMR
nlast = nthermr2
commands = commands.format(**locals())
# acer
nthermal_acer_in = nlast
for i, temperature in enumerate(temperatures):
# Extend input with an ACER run for each temperature
nace = nthermal_acer_in + 1 + 2*i
ndir = nace + 1
ext = '{:02}'.format(i + 1)
commands += _THERMAL_TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = output_dir / "ace_{:.1f}".format(temperature)
tapeout[ndir] = output_dir / "xsdir_{:.1f}".format(temperature)
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
ace = output_dir / ace
xsdir = (ace.parent / "xsdir") if xsdir is None else Path(xsdir)
with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
# Concatenate ACE and xsdir files together
for temperature in temperatures:
ace_in = output_dir / "ace_{:.1f}".format(temperature)
ace_file.write(ace_in.read_text())
xsdir_in = output_dir / "xsdir_{:.1f}".format(temperature)
xsdir_file.write(xsdir_in.read_text())
# Remove ACE/xsdir files for each temperature
for temperature in temperatures:
(output_dir / "ace_{:.1f}".format(temperature)).unlink()
(output_dir / "xsdir_{:.1f}".format(temperature)).unlink()