diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 9c73152d45..9024b932e4 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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: diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index c8d630851b..7993ec4004 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -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, diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 05b6b8d2d6..696310ccaf 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -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)}) diff --git a/src/reaction.cpp b/src/reaction.cpp index b733b20283..31bf8efd33 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -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) { diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 3cf036aa2c..40a286331e 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -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)