Skip to content

Commit

Permalink
Correct normalization of thermal elastic in non standard ENDF-6 files (…
Browse files Browse the repository at this point in the history 
…#3234)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
  • Loading branch information
@marquezj @paulromano
marquezj and paulromano authored Feb 12, 2025
1 parent 0439320 commit e9ddf88
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 58 deletions.
137 changes: 82 additions & 55 deletions openmc/data/thermal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
from io import StringIO
import itertools
import os
import re
import tempfile
from warnings import warn

Expand Down Expand Up @@ -871,7 +870,7 @@ def from_ace(cls, ace_or_filename, name=None):
@classmethod
def from_njoy(cls, filename, filename_thermal, temperatures=None,
evaluation=None, evaluation_thermal=None,
use_endf_data=True, **kwargs):
use_endf_data=True, divide_incoherent_elastic=False, **kwargs):
"""Generate thermal scattering data by running NJOY.
Parameters
Expand All @@ -894,8 +893,13 @@ def from_njoy(cls, filename, filename_thermal, temperatures=None,
use_endf_data : bool
If the material has incoherent elastic scattering, the ENDF data
will be used rather than the ACE data.
divide_incoherent_elastic : bool
Divide incoherent elastic cross section by number of principal
atoms. This is not part of the ENDF-6 standard but it is how it is
processed by NJOY.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal`
Keyword arguments passed to
:func:`openmc.data.njoy.make_ace_thermal`
Returns
-------
Expand All @@ -920,7 +924,7 @@ def from_njoy(cls, filename, filename_thermal, temperatures=None,

# Load ENDF data to replace incoherent elastic
if use_endf_data:
data_endf = cls.from_endf(filename_thermal)
data_endf = cls.from_endf(filename_thermal, divide_incoherent_elastic)
if data_endf.elastic is not None:
# Get appropriate temperatures
if temperatures is None:
Expand All @@ -938,14 +942,18 @@ def from_njoy(cls, filename, filename_thermal, temperatures=None,
return data

@classmethod
def from_endf(cls, ev_or_filename):
def from_endf(cls, ev_or_filename, divide_incoherent_elastic=False):
"""Generate thermal scattering data from an ENDF file
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
divide_incoherent_elastic : bool
Divide incoherent elastic cross section by number of principal
atoms. This is not part of the ENDF-6 standard but it is how it is
processed by NJOY.
Returns
-------
Expand All @@ -958,56 +966,6 @@ def from_endf(cls, ev_or_filename):
else:
ev = endf.Evaluation(ev_or_filename)

# Read coherent/incoherent elastic data
elastic = None
if (7, 2) in ev.section:
# Define helper functions to avoid duplication
def get_coherent_elastic(file_obj):
# Get structure factor at first temperature
params, S = endf.get_tab1_record(file_obj)
strT = _temperature_str(params[0])
n_temps = params[2]
bragg_edges = S.x
xs = {strT: CoherentElastic(bragg_edges, S.y)}
distribution = {strT: CoherentElasticAE(xs[strT])}

# Get structure factor for subsequent temperatures
for _ in range(n_temps):
params, S = endf.get_list_record(file_obj)
strT = _temperature_str(params[0])
xs[strT] = CoherentElastic(bragg_edges, S)
distribution[strT] = CoherentElasticAE(xs[strT])
return xs, distribution

def get_incoherent_elastic(file_obj):
params, W = endf.get_tab1_record(file_obj)
bound_xs = params[0]
xs = {}
distribution = {}
for T, debye_waller in zip(W.x, W.y):
strT = _temperature_str(T)
xs[strT] = IncoherentElastic(bound_xs, debye_waller)
distribution[strT] = IncoherentElasticAE(debye_waller)
return xs, distribution

file_obj = StringIO(ev.section[7, 2])
lhtr = endf.get_head_record(file_obj)[2]
if lhtr == 1:
# coherent elastic
xs, distribution = get_coherent_elastic(file_obj)
elif lhtr == 2:
# incoherent elastic
xs, distribution = get_incoherent_elastic(file_obj)
elif lhtr == 3:
# mixed coherent / incoherent elastic
xs_c, dist_c = get_coherent_elastic(file_obj)
xs_i, dist_i = get_incoherent_elastic(file_obj)
assert sorted(xs_c) == sorted(xs_i)
xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c}
distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c}

elastic = ThermalScatteringReaction(xs, distribution)

# Read incoherent inelastic data
assert (7, 4) in ev.section, 'No MF=7, MT=4 found in thermal scattering'
file_obj = StringIO(ev.section[7, 4])
Expand All @@ -1022,6 +980,7 @@ def get_incoherent_elastic(file_obj):
data['A0'] = awr = B[2]
data['e_max'] = energy_max = B[3]
data['M0'] = B[5]
free_xs = data['free_atom_xs'] / data['M0']

# Get information about non-principal atoms
n_non_principal = params[5]
Expand Down Expand Up @@ -1066,6 +1025,74 @@ def get_incoherent_elastic(file_obj):
_, Teff = endf.get_tab1_record(file_obj)
data['effective_temperature'].append(Teff)

# Read coherent/incoherent elastic data
elastic = None
if (7, 2) in ev.section:
# Define helper functions to avoid duplication
def get_coherent_elastic(file_obj):
# Get structure factor at first temperature
params, S = endf.get_tab1_record(file_obj)
strT = _temperature_str(params[0])
n_temps = params[2]
bragg_edges = S.x
xs = {strT: CoherentElastic(bragg_edges, S.y)}
distribution = {strT: CoherentElasticAE(xs[strT])}

# Get structure factor for subsequent temperatures
for _ in range(n_temps):
params, S = endf.get_list_record(file_obj)
strT = _temperature_str(params[0])
xs[strT] = CoherentElastic(bragg_edges, S)
distribution[strT] = CoherentElasticAE(xs[strT])
return xs, distribution

def get_incoherent_elastic(file_obj, natom):
params, W = endf.get_tab1_record(file_obj)
bound_xs = params[0]/natom

# Check whether divide_incoherent_elastic was applied correctly
if abs(free_xs - bound_xs/(1 + 1/data['A0'])**2) > 0.5:
if divide_incoherent_elastic:
msg = (
'Thermal scattering evaluation follows ENDF-6 '
'definition of bound cross section but '
'divide_incoherent_elastic=True.'
)
else:
msg = (
'Thermal scattering evaluation follows NJOY '
'definition of bound cross section but '
'divide_incoherent_elastic=False.'
)
warn(msg)

xs = {}
distribution = {}
for T, debye_waller in zip(W.x, W.y):
strT = _temperature_str(T)
xs[strT] = IncoherentElastic(bound_xs, debye_waller)
distribution[strT] = IncoherentElasticAE(debye_waller)
return xs, distribution

file_obj = StringIO(ev.section[7, 2])
lhtr = endf.get_head_record(file_obj)[2]
natom = data['M0'] if divide_incoherent_elastic else 1
if lhtr == 1:
# coherent elastic
xs, distribution = get_coherent_elastic(file_obj)
elif lhtr == 2:
# incoherent elastic
xs, distribution = get_incoherent_elastic(file_obj, natom)
elif lhtr == 3:
# mixed coherent / incoherent elastic
xs_c, dist_c = get_coherent_elastic(file_obj)
xs_i, dist_i = get_incoherent_elastic(file_obj, natom)
assert sorted(xs_c) == sorted(xs_i)
xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c}
distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c}

elastic = ThermalScatteringReaction(xs, distribution)

name = ev.target['zsymam'].strip()
instance = cls(name, awr, energy_max, kTs)
if elastic is not None:
Expand Down
6 changes: 3 additions & 3 deletions tests/unit_tests/test_data_thermal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ def hzrh():
"""H in ZrH thermal scattering data."""
endf_data = os.environ['OPENMC_ENDF_DATA']
filename = os.path.join(endf_data, 'thermal_scatt', 'tsl-HinZrH.endf')
return openmc.data.ThermalScattering.from_endf(filename)
return openmc.data.ThermalScattering.from_endf(filename, divide_incoherent_elastic=True)


@pytest.fixture(scope='module')
Expand All @@ -64,7 +64,7 @@ def sio2():
"""SiO2 thermal scattering data."""
endf_data = os.environ['OPENMC_ENDF_DATA']
filename = os.path.join(endf_data, 'thermal_scatt', 'tsl-SiO2.endf')
return openmc.data.ThermalScattering.from_endf(filename)
return openmc.data.ThermalScattering.from_endf(filename, divide_incoherent_elastic=True)


def test_h2o_attributes(h2o):
Expand Down Expand Up @@ -144,7 +144,7 @@ def test_continuous_dist(h2o_njoy):
def test_h2o_endf():
endf_data = os.environ['OPENMC_ENDF_DATA']
filename = os.path.join(endf_data, 'thermal_scatt', 'tsl-HinH2O.endf')
h2o = openmc.data.ThermalScattering.from_endf(filename)
h2o = openmc.data.ThermalScattering.from_endf(filename, divide_incoherent_elastic=True)
assert not h2o.elastic
assert h2o.atomic_weight_ratio == pytest.approx(0.99917)
assert h2o.energy_max == pytest.approx(3.99993)
Expand Down

0 comments on commit e9ddf88

Please sign in to comment.