From d81aaeca981e2e10d06adfc560b6726d979f9a6e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 9 Aug 2019 11:25:16 -0500 Subject: [PATCH 1/9] Store fission heating IncidentNeutron.from_ace IncidentNeutron.from_ace now computes the fission heating and a fission-less heating coefficient. The fission heating is the product of the heating number and the fission cross section, and stored as MT318. The fission-less heating is the heating from all reactions except fission, computed as heating_number * (total_xs - fission_xs). The MT number for this is taken to be 999 as a temporary value. A test is added for Am244 in test_data_neutron.py to examine the new heating values --- openmc/data/neutron.py | 12 ++++++++++++ tests/unit_tests/test_data_neutron.py | 12 ++++++++++++ 2 files changed, 24 insertions(+) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 9c73152d45..e903b6a92f 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -633,6 +633,18 @@ class IncidentNeutron(EqualityMixin): rx = Reaction.from_ace(ace, i) data.reactions[rx.mt] = rx + # If present, use fission xs to compute "fission heating" coefficient + fission_reaction = data.reactions.get(18) + if fission_reaction is not None: + fission_xs = fission_reaction.xs[strT].y + + # Compute "fission-less" heating coefficient + no_fission_heating = Reaction(999) + no_fission_heating.xs[strT] = Tabulated1D( + energy, heating_number * (total_xs - fission_xs)) + no_fission_heating.redundant = True + data.reactions[999] = no_fission_heating + # Some photon production reactions may be assigned to MTs that don't # exist, usually MT=4. In this case, we create a new reaction and add # them diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 3cf036aa2c..7fccb33d88 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -205,6 +205,18 @@ def test_derived_products(am244): assert total_neutron.yield_(6e6) == pytest.approx(4.2558) +def test_heating(run_in_tmpdir, am244): + assert 999 in am244.reactions # TBD + strT = min(am244.reactions[1].xs.keys()) + total_xs = am244.reactions[1].xs[strT].y + fission_xs = am244.reactions[18].xs[strT].y + total_heating = am244.reactions[301].xs[strT].y + no_fission_heating = am244.reactions[999].xs[strT].y + heating_number = total_heating / total_xs + assert no_fission_heating == pytest.approx( + heating_number * (total_xs - fission_xs)) + + def test_urr(pu239): for T, ptable in pu239.urr.items(): assert T.endswith('K') From 81c6259289da9f661096e6b91a7bec9460a30a3b Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 19 Aug 2019 11:39:11 -0500 Subject: [PATCH 2/9] Write neutron fission-less heating coefficient to HDF5 Mark reaction 999 as one to be kept during the export process --- openmc/data/neutron.py | 2 +- tests/unit_tests/test_data_neutron.py | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index e903b6a92f..a53194aa78 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -453,7 +453,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, 444, 999) if not (photon_rx or rx.mt in keep_mts): continue diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 7fccb33d88..22aebe20ad 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -206,7 +206,7 @@ def test_derived_products(am244): def test_heating(run_in_tmpdir, am244): - assert 999 in am244.reactions # TBD + assert 999 in am244.reactions strT = min(am244.reactions[1].xs.keys()) total_xs = am244.reactions[1].xs[strT].y fission_xs = am244.reactions[18].xs[strT].y @@ -215,6 +215,12 @@ def test_heating(run_in_tmpdir, am244): heating_number = total_heating / total_xs assert no_fission_heating == pytest.approx( heating_number * (total_xs - fission_xs)) + # Re-read to ensure data is written + am244.export_to_hdf5("Am244.h5") + new = openmc.data.IncidentNeutron.from_hdf5("Am244.h5") + assert 999 in new.reactions + new_rxn = new.reactions[999].xs[strT].y + assert (new_rxn == no_fission_heating).all() def test_urr(pu239): From 7badc48c225671a087fa78300bfd3dbcd16f2906 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 3 Sep 2019 16:15:14 -0500 Subject: [PATCH 3/9] Instruct njoy interface to retain heatr output HEATR contains the 318 fission heating data needed to construct the energy deposition tally. Until now, the output tape was left in the njoy temporary directory with the other output tapes. This commit instructs openmc.data.njoy to move the heatr output from the running directory used in njoy.run to the working directory used in njoy.make_ace Removed unused imports argparse and sys --- openmc/data/neutron.py | 5 ++--- openmc/data/njoy.py | 8 ++++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index a53194aa78..22c9d03054 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -813,9 +813,8 @@ 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')) + for key in ["ace", "xsdir", "pendf", "heatr"]: + kwargs.setdefault(key, os.path.join(tmpdir, key)) kwargs['evaluation'] = evaluation make_ace(filename, temperatures, **kwargs) diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index c8d630851b..3fb75751f5 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -1,10 +1,8 @@ -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 . import endf @@ -237,8 +235,9 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, 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 + 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, optional Indicating whether to add gas production data when running NJOY purr : bool, optional @@ -292,6 +291,7 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, if heatr: nheatr_in = nlast nheatr = nheatr_in + 1 + tapeout[nheatr] = "heatr" if heatr is True else heatr commands += _TEMPLATE_HEATR nlast = nheatr From 2675ebc2c6bb22f36c791280dc6332d8cb9227e2 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 3 Sep 2019 17:09:15 -0500 Subject: [PATCH 4/9] Compute fission heating and fission-less heating with njoy This reverts d81aaeca9 and 81c625928 in the following ways. Fission heating data, MT318, is pulled from the heatr file produced when running NJOY. Reaction data that is potentially temperature dependent is set onto the IncidentNeutron object by scaling MT318 by the ratio of the fission cross section used in HEATR and other fission cross sections stored on the object. This produces potentially many MT318 fission heating KERMA coefficients on the nuclide. The fission-less heating coefficient, MT999, if computed by subtracting MT318 from MT301, total heating coefficients. This reaction is marked as redundant, as it can easily be reconstructed. MT318 is allowed to be written to the HDF5 file, while 999 is not. When reading back in the library, MT999 is rebuilt in exactly the same manner. The test_heating test in tests/unit_tests/test_data_neutron.py has been updated given the changes in this commit. It is worth noting that the heating coefficients computed in NJOY only contain prompt neutrons, as k_{i,j}(E) = sigma_{i,j}(E) * (E + Q - \bar{E}) where k_{i,j} is the kerma coefficient for reaction j of material i, sigma_{i,j} is the corresponding reaction cross section, E is the energy of incident particle, Q is the mass-difference Q value, and \bar{E} is the average energy of secondary particles. Source: NJOY16 Manual on HEATR --- openmc/data/neutron.py | 71 +++++++++++++++++++++------ tests/unit_tests/test_data_neutron.py | 26 +++++----- 2 files changed, 67 insertions(+), 30 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 22c9d03054..ae3b2c6c22 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 @@ -453,7 +451,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, 999) + 203, 204, 205, 206, 207, 301, 318, 444) if not (photon_rx or rx.mt in keep_mts): continue @@ -554,6 +552,22 @@ class IncidentNeutron(EqualityMixin): fer_group = group['fission_energy_release'] data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) + # Rebuild fission-less 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: + fission_less_heating = Reaction(999) + fission_less_heating.redundant = True + for strT, total in total_heating.xs.items(): + fission = fission_heating.xs.get(strT) + if fission is None: + continue + fission_less_heating.xs[strT] = Tabulated1D( + total.x, total.y - fission(total.x), + breakpoints=total.breakpoints, + interpolation=total.interpolation) + data.reactions[999] = fission_less_heating + return data @classmethod @@ -633,18 +647,6 @@ class IncidentNeutron(EqualityMixin): rx = Reaction.from_ace(ace, i) data.reactions[rx.mt] = rx - # If present, use fission xs to compute "fission heating" coefficient - fission_reaction = data.reactions.get(18) - if fission_reaction is not None: - fission_xs = fission_reaction.xs[strT].y - - # Compute "fission-less" heating coefficient - no_fission_heating = Reaction(999) - no_fission_heating.xs[strT] = Tabulated1D( - energy, heating_number * (total_xs - fission_xs)) - no_fission_heating.redundant = True - data.reactions[999] = no_fission_heating - # Some photon production reactions may be assigned to MTs that don't # exist, usually MT=4. In this case, we create a new reaction and add # them @@ -828,6 +830,43 @@ 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 + heatr = Evaluation(kwargs["heatr"]) + f318 = StringIO(heatr.section[3, 318]) + get_head_record(f318) + _params, fission_kerma = get_tab1_record(f318) + # Assume that only cross section changes with temperature + # to compute heating data for multiple temperatures + f18 = StringIO(heatr.section[3, 18]) + get_head_record(f18) + _params, fission_xs_heatr = get_tab1_record(f18) + heat_num = Tabulated1D( + fission_kerma.x, + fission_kerma.y / fission_xs_heatr(fission_kerma.x), + breakpoints=fission_kerma.breakpoints, + interpolation=fission_kerma.interpolation) + + fission_heating = Reaction(318) + fission_less_heating = Reaction(999) + fission_less_heating.redundant = True + + for strT, fission_xs in data.reactions[18].xs.items(): + total_heating = data.reactions[301].xs.get(strT) + if total_heating is None: + continue + heater_at_temp = Tabulated1D( + fission_xs.x, heat_num(fission_xs.x) * fission_xs.y, + breakpoints=fission_xs.breakpoints, + interpolation=fission_xs.interpolation) + fission_heating.xs[strT] = heater_at_temp + fission_less_heating.xs[strT] = Tabulated1D( + total_heating.x, + total_heating.y - heater_at_temp(total_heating.x), + breakpoints=total_heating.breakpoints, + interpolation=total_heating.interpolation) + + data.reactions[318] = fission_heating + data.reactions[999] = fission_less_heating # Add 0K elastic scattering cross section if '0K' not in data.energy: diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 22aebe20ad..c5da4acfaf 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -206,21 +206,19 @@ def test_derived_products(am244): def test_heating(run_in_tmpdir, am244): + assert 318 in am244.reactions assert 999 in am244.reactions - strT = min(am244.reactions[1].xs.keys()) - total_xs = am244.reactions[1].xs[strT].y - fission_xs = am244.reactions[18].xs[strT].y - total_heating = am244.reactions[301].xs[strT].y - no_fission_heating = am244.reactions[999].xs[strT].y - heating_number = total_heating / total_xs - assert no_fission_heating == pytest.approx( - heating_number * (total_xs - fission_xs)) - # Re-read to ensure data is written - am244.export_to_hdf5("Am244.h5") - new = openmc.data.IncidentNeutron.from_hdf5("Am244.h5") - assert 999 in new.reactions - new_rxn = new.reactions[999].xs[strT].y - assert (new_rxn == no_fission_heating).all() + # 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): From a0bac56ebdb5a7062a8d6ca0f115b1ee0a195420 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Sep 2019 13:47:32 -0500 Subject: [PATCH 5/9] Add fission-heating and non-fission-heating to reactions MT318 and MT999 for fission heating and non-fission heating have been added to the reactions supported in openmc/data/reaction.py and src/reaction.cpp Changed internal language from fission-less heating to non-fission heating --- openmc/data/neutron.py | 16 ++++++++-------- openmc/data/reaction.py | 1 + src/reaction.cpp | 4 ++++ 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index ae3b2c6c22..0f474446cf 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -556,17 +556,17 @@ class IncidentNeutron(EqualityMixin): total_heating = data.reactions.get(301) fission_heating = data.reactions.get(318) if total_heating is not None and fission_heating is not None: - fission_less_heating = Reaction(999) - fission_less_heating.redundant = True + 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 - fission_less_heating.xs[strT] = Tabulated1D( + non_fission_heating.xs[strT] = Tabulated1D( total.x, total.y - fission(total.x), breakpoints=total.breakpoints, interpolation=total.interpolation) - data.reactions[999] = fission_less_heating + data.reactions[999] = non_fission_heating return data @@ -847,8 +847,8 @@ class IncidentNeutron(EqualityMixin): interpolation=fission_kerma.interpolation) fission_heating = Reaction(318) - fission_less_heating = Reaction(999) - fission_less_heating.redundant = True + non_fission_heating = Reaction(999) + non_fission_heating.redundant = True for strT, fission_xs in data.reactions[18].xs.items(): total_heating = data.reactions[301].xs.get(strT) @@ -859,14 +859,14 @@ class IncidentNeutron(EqualityMixin): breakpoints=fission_xs.breakpoints, interpolation=fission_xs.interpolation) fission_heating.xs[strT] = heater_at_temp - fission_less_heating.xs[strT] = Tabulated1D( + non_fission_heating.xs[strT] = Tabulated1D( total_heating.x, total_heating.y - heater_at_temp(total_heating.x), breakpoints=total_heating.breakpoints, interpolation=total_heating.interpolation) data.reactions[318] = fission_heating - data.reactions[999] = fission_less_heating + data.reactions[999] = non_fission_heating # Add 0K elastic scattering cross section if '0K' not in data.energy: 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) { From 6b23c21327ea660f5bafabc4e192e1b72252bbda Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Sep 2019 13:28:06 -0500 Subject: [PATCH 6/9] Pass output_dir to make_ace to collect files The user can now pass a single argument indicating where output tapes for each njoy module should be written. The default argument is the current working directory, useful for removing any temporary files created during testing. The pendf, broadr, heatr, gaspr, and purr arguments can either be booleans or strings. A boolean value indicates the module output should be written into a file named after the module, e.g. heatr output written to /heatr. ace and xsdir can be None or strings. A value of None means the final ace and xsdir files will be written to /ace and /xsdir, respectively. The temporary temperature-dependent ace and xsdir files are now always named "ace_" and "xsdir_", as the ace and xsdir arguments can be more than "ace" and "xsdir". openmc.data.IncidentNeutron.from_njoy knows to use this output_dir argument to ensure all created tape files are written into the temporary directory or the user requested directory, if applicable. This commit was written because simply writing the broadr, gaspr, and purr output files to "broadr", "gaspr", and "purr" resulted in left over files after running tests. --- openmc/data/neutron.py | 6 ++-- openmc/data/njoy.py | 76 ++++++++++++++++++++++++++++++------------ 2 files changed, 59 insertions(+), 23 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 0f474446cf..26c0a71c5a 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -815,8 +815,10 @@ class IncidentNeutron(EqualityMixin): """ with tempfile.TemporaryDirectory() as tmpdir: # Run NJOY to create an ACE library - for key in ["ace", "xsdir", "pendf", "heatr"]: - kwargs.setdefault(key, os.path.join(tmpdir, key)) + kwargs.setdefault("output_dir", tmpdir) + for key in ("ace", "xsdir", "pendf", "heatr", "broadr", + "gaspr", "purr"): + kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key)) kwargs['evaluation'] = evaluation make_ace(filename, temperatures, **kwargs) diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 3fb75751f5..a648a53515 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -213,11 +213,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, ace=None, xsdir=None, + output_dir=None, pendf=False, 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 + File names can be passed to + ``[ace, 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 @@ -226,22 +236,30 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, 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 + 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"`` + 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 + 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, optional - Indicating whether to add gas production data when running NJOY - purr : bool, optional - Indicating whether to add probability table when running NJOY + 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. acer : bool, optional Indicating whether to generate ACE file when running NJOY evaluation : openmc.data.endf.Evaluation, optional @@ -256,6 +274,12 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, If the NJOY process returns with a non-zero status """ + if output_dir is None: + output_dir = os.path.abspath(os.curdir) + else: + assert os.path.isdir(output_dir), output_dir + output_dir = os.path.abspath(output_dir) + ev = evaluation if evaluation is not None else endf.Evaluation(filename) mat = ev.material zsymam = ev.target['zsymam'] @@ -274,8 +298,9 @@ 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] = (os.path.join(output_dir, "pendf") if pendf is True + else pendf) # reconr commands += _TEMPLATE_RECONR @@ -284,6 +309,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, # broadr if broadr: nbroadr = nlast + 1 + tapeout[nbroadr] = ( + os.path.join(output_dir, "broadr") if broadr is True else broadr) commands += _TEMPLATE_BROADR nlast = nbroadr @@ -291,7 +318,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, if heatr: nheatr_in = nlast nheatr = nheatr_in + 1 - tapeout[nheatr] = "heatr" if heatr is True else heatr + tapeout[nheatr] = (os.path.join(output_dir, "heatr") if heatr is True + else heatr) commands += _TEMPLATE_HEATR nlast = nheatr @@ -299,6 +327,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, if gaspr: ngaspr_in = nlast ngaspr = ngaspr_in + 1 + tapeout[ngaspr] = (os.path.join(output_dir, "gaspr") if gaspr is True + else gaspr) commands += _TEMPLATE_GASPR nlast = ngaspr @@ -306,6 +336,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, if purr: npurr_in = nlast npurr = npurr_in + 1 + tapeout[npurr] = (os.path.join(output_dir, "purr") if purr is True + else purr) commands += _TEMPLATE_PURR nlast = npurr @@ -323,16 +355,18 @@ 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: + ace = os.path.join(output_dir, "ace") if ace is None else ace + xsdir = os.path.join(output_dir, "xsdir") if xsdir is None else xsdir with open(ace, 'w') as ace_file, open(xsdir, '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 @@ -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, From c7a8b6606e444046cd8b04609c74a549650ff18e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Sep 2019 14:46:32 -0500 Subject: [PATCH 7/9] Better collection of fission heating at temperatures Use openmc.data.endf.get_evaluations to read in potentially many 318 evaluations from the heatr file. The evaluated temperatures are pulled using the target attribute, rounded to the nearest integer for consistency, e.g. 293.6K -> 294K. --- openmc/data/neutron.py | 47 +++++++++++++++--------------------------- 1 file changed, 17 insertions(+), 30 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 26c0a71c5a..37230e62ca 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -13,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 @@ -833,39 +834,25 @@ class IncidentNeutron(EqualityMixin): if (1, 458) in ev.section: data.fission_energy = FissionEnergyRelease.from_endf(ev, data) # Add 318 fission heating data from heatr - heatr = Evaluation(kwargs["heatr"]) - f318 = StringIO(heatr.section[3, 318]) - get_head_record(f318) - _params, fission_kerma = get_tab1_record(f318) - # Assume that only cross section changes with temperature - # to compute heating data for multiple temperatures - f18 = StringIO(heatr.section[3, 18]) - get_head_record(f18) - _params, fission_xs_heatr = get_tab1_record(f18) - heat_num = Tabulated1D( - fission_kerma.x, - fission_kerma.y / fission_xs_heatr(fission_kerma.x), - breakpoints=fission_kerma.breakpoints, - interpolation=fission_kerma.interpolation) - - fission_heating = Reaction(318) non_fission_heating = Reaction(999) non_fission_heating.redundant = True + fission_heating = Reaction(318) - for strT, fission_xs in data.reactions[18].xs.items(): - total_heating = data.reactions[301].xs.get(strT) - if total_heating is None: + 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 - heater_at_temp = Tabulated1D( - fission_xs.x, heat_num(fission_xs.x) * fission_xs.y, - breakpoints=fission_xs.breakpoints, - interpolation=fission_xs.interpolation) - fission_heating.xs[strT] = heater_at_temp - non_fission_heating.xs[strT] = Tabulated1D( - total_heating.x, - total_heating.y - heater_at_temp(total_heating.x), - breakpoints=total_heating.breakpoints, - interpolation=total_heating.interpolation) + 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 From 05f33849f51f4d6fba1d75979d958e1a24a9d0f8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 5 Sep 2019 08:35:52 -0500 Subject: [PATCH 8/9] Remove breakpoints, interpolation for non-fission heating NJOY handles linearization of the cross sections in RECONR: Section 3.2 NJOY2016 Manual, LA-UR-17-20093 --- openmc/data/neutron.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 37230e62ca..16ff032a6d 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -553,7 +553,7 @@ class IncidentNeutron(EqualityMixin): fer_group = group['fission_energy_release'] data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) - # Rebuild fission-less heating + # 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: @@ -564,9 +564,7 @@ class IncidentNeutron(EqualityMixin): if fission is None: continue non_fission_heating.xs[strT] = Tabulated1D( - total.x, total.y - fission(total.x), - breakpoints=total.breakpoints, - interpolation=total.interpolation) + total.x, total.y - fission(total.x)) data.reactions[999] = non_fission_heating return data From 76adfb231c53ac276643046c8a5fbd728ea194f9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 5 Sep 2019 09:18:28 -0500 Subject: [PATCH 9/9] Combine acer, ace arguments to openmc.data.make_ace Passing ace to make_ace is not supported. Use acer instead. The acer argument now signals if the acer njoy module should be run and where to save the output file. xsdir defaults to None now, rather than "xsdir". This change allows the xsdir file to be written to the same directory as acer, or a user supplied location potentially separate from ace file. --- openmc/data/neutron.py | 5 ++- openmc/data/njoy.py | 52 +++++++++++++-------------- tests/unit_tests/test_data_neutron.py | 2 +- 3 files changed, 29 insertions(+), 30 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 16ff032a6d..9024b932e4 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -815,14 +815,13 @@ class IncidentNeutron(EqualityMixin): with tempfile.TemporaryDirectory() as tmpdir: # Run NJOY to create an ACE library kwargs.setdefault("output_dir", tmpdir) - for key in ("ace", "xsdir", "pendf", "heatr", "broadr", - "gaspr", "purr"): + 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) diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index a648a53515..7993ec4004 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -4,6 +4,7 @@ import os import shutil from subprocess import Popen, PIPE, STDOUT, CalledProcessError import tempfile +from pathlib import Path from . import endf @@ -213,14 +214,14 @@ 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=None, xsdir=None, +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, acer=True, evaluation=None, + heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs): """Generate incident neutron ACE file from an ENDF file File names can be passed to - ``[ace, xsdir, pendf, broadr, heatr, gaspr, purr]`` + ``[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 @@ -235,10 +236,13 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=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. Defaults to ``"ace"`` + 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"`` + 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]`` @@ -260,8 +264,6 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, purr : bool or str, optional Indicating whether to add probability table when running NJOY. If string, write the output tape to this file. - 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. @@ -272,13 +274,16 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=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 = os.path.abspath(os.curdir) + output_dir = Path() else: - assert os.path.isdir(output_dir), output_dir - output_dir = os.path.abspath(output_dir) + 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 @@ -299,8 +304,7 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, tapein = {nendf: filename} tapeout = {} if pendf: - tapeout[npendf] = (os.path.join(output_dir, "pendf") if pendf is True - else pendf) + tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf # reconr commands += _TEMPLATE_RECONR @@ -309,8 +313,7 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, # broadr if broadr: nbroadr = nlast + 1 - tapeout[nbroadr] = ( - os.path.join(output_dir, "broadr") if broadr is True else broadr) + tapeout[nbroadr] = (output_dir / "broadr") if broadr is True else broadr commands += _TEMPLATE_BROADR nlast = nbroadr @@ -318,8 +321,7 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, if heatr: nheatr_in = nlast nheatr = nheatr_in + 1 - tapeout[nheatr] = (os.path.join(output_dir, "heatr") if heatr is True - else heatr) + tapeout[nheatr] = (output_dir / "heatr") if heatr is True else heatr commands += _TEMPLATE_HEATR nlast = nheatr @@ -327,8 +329,7 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, if gaspr: ngaspr_in = nlast ngaspr = ngaspr_in + 1 - tapeout[ngaspr] = (os.path.join(output_dir, "gaspr") if gaspr is True - else gaspr) + tapeout[ngaspr] = (output_dir / "gaspr") if gaspr is True else gaspr commands += _TEMPLATE_GASPR nlast = ngaspr @@ -336,8 +337,7 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, if purr: npurr_in = nlast npurr = npurr_in + 1 - tapeout[npurr] = (os.path.join(output_dir, "purr") if purr is True - else purr) + tapeout[npurr] = (output_dir / "purr") if purr is True else purr commands += _TEMPLATE_PURR nlast = npurr @@ -361,15 +361,15 @@ def make_ace(filename, temperatures=None, ace=None, xsdir=None, run(commands, tapein, tapeout, **kwargs) if acer: - ace = os.path.join(output_dir, "ace") if ace is None else ace - xsdir = os.path.join(output_dir, "xsdir") if xsdir is None else xsdir - 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() - # 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: diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index c5da4acfaf..40a286331e 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -481,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)