mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Correct normalization of thermal elastic in non standard ENDF-6 files (#3234)
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Jose Ignacio Marquez Damian
GitHub
parent
04393200cc
commit
e9ddf885e7
2 changed files with 85 additions and 58 deletions
|
|
@ -5,7 +5,6 @@ from numbers import Real
|
|||
from io import StringIO
|
||||
import itertools
|
||||
import os
|
||||
import re
|
||||
import tempfile
|
||||
from warnings import warn
|
||||
|
||||
|
|
@ -871,7 +870,7 @@ class ThermalScattering(EqualityMixin):
|
|||
@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
|
||||
|
|
@ -894,8 +893,13 @@ class ThermalScattering(EqualityMixin):
|
|||
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
|
||||
-------
|
||||
|
|
@ -920,7 +924,7 @@ class ThermalScattering(EqualityMixin):
|
|||
|
||||
# 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:
|
||||
|
|
@ -938,7 +942,7 @@ class ThermalScattering(EqualityMixin):
|
|||
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
|
||||
|
|
@ -946,6 +950,10 @@ class ThermalScattering(EqualityMixin):
|
|||
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
|
||||
-------
|
||||
|
|
@ -958,56 +966,6 @@ class ThermalScattering(EqualityMixin):
|
|||
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])
|
||||
|
|
@ -1022,6 +980,7 @@ class ThermalScattering(EqualityMixin):
|
|||
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]
|
||||
|
|
@ -1066,6 +1025,74 @@ class ThermalScattering(EqualityMixin):
|
|||
_, 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:
|
||||
|
|
|
|||
|
|
@ -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')
|
||||
|
|
@ -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):
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue