From 2675ebc2c6bb22f36c791280dc6332d8cb9227e2 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 3 Sep 2019 17:09:15 -0500 Subject: [PATCH] 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):