mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #2134 from paulromano/thermal-mixed-elastic
Add support for mixed incoherent/coherent elastic thermal scattering
This commit is contained in:
Amelia J Trainer
GitHub
commit
a18abd4fd9
15 changed files with 476 additions and 48 deletions
|
|
@ -339,6 +339,16 @@ Incoherent elastic scattering
|
|||
[eV\ :math:`^{-1}`].
|
||||
:Attributes: - **type** (*char[]*) -- 'IncoherentElastic'
|
||||
|
||||
Sum of functions
|
||||
----------------
|
||||
|
||||
:Object type: Group
|
||||
:Attributes: - **type** (*char[]*) -- "Sum"
|
||||
- **n** (*int*) -- Number of functions
|
||||
:Datasets:
|
||||
- ***func_<i>** (:ref:`function <1d_functions>`) -- Dataset for the
|
||||
i-th function (indexing starts at 1)
|
||||
|
||||
.. _angle_energy:
|
||||
|
||||
--------------------------
|
||||
|
|
@ -501,6 +511,19 @@ equiprobable bins.
|
|||
- **skewed** (*int8_t*) -- Whether discrete angles are equi-probable
|
||||
(0) or have a skewed distribution (1).
|
||||
|
||||
Mixed Elastic
|
||||
-------------
|
||||
|
||||
This angle-energy distribution is used when an evaluation specifies both
|
||||
coherent and incoherent elastic thermal neutron scattering.
|
||||
|
||||
:Object type: Group
|
||||
:Attributes: - **type** (*char[]*) -- "mixed_elastic"
|
||||
:Groups: - **coherent** -- Distribution for coherent elastic scattering. The
|
||||
format is given in :ref:`angle_energy`.
|
||||
- **incoherent** -- Distribution for incoherent elastic scattering.
|
||||
The format is given in :ref:`angle_energy`.
|
||||
|
||||
.. _energy_distribution:
|
||||
|
||||
--------------------
|
||||
|
|
|
|||
|
|
@ -116,6 +116,7 @@ Angle-Energy Distributions
|
|||
IncoherentElasticAE
|
||||
IncoherentElasticAEDiscrete
|
||||
IncoherentInelasticAEDiscrete
|
||||
MixedElasticAE
|
||||
|
||||
Resonance Data
|
||||
--------------
|
||||
|
|
|
|||
|
|
@ -126,6 +126,26 @@ private:
|
|||
debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1]
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Sum of multiple 1D functions
|
||||
//==============================================================================
|
||||
|
||||
class Sum1D : public Function1D {
|
||||
public:
|
||||
// Constructors
|
||||
explicit Sum1D(hid_t group);
|
||||
|
||||
//! Evaluate each function and sum results
|
||||
//! \param[in] x independent variable
|
||||
//! \return Function evaluated at x
|
||||
double operator()(double E) const override;
|
||||
|
||||
const unique_ptr<Function1D>& functions(int i) const { return functions_[i]; }
|
||||
|
||||
private:
|
||||
vector<unique_ptr<Function1D>> functions_; //!< individual functions
|
||||
};
|
||||
|
||||
//! Read 1D function from HDF5 dataset
|
||||
//! \param[in] group HDF5 group containing dataset
|
||||
//! \param[in] name Name of dataset
|
||||
|
|
|
|||
|
|
@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false);
|
|||
vector<std::string> group_names(hid_t group_id);
|
||||
vector<hsize_t> object_shape(hid_t obj_id);
|
||||
std::string object_name(hid_t obj_id);
|
||||
hid_t open_object(hid_t group_id, const std::string& name);
|
||||
void close_object(hid_t obj_id);
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility functions
|
||||
|
|
|
|||
|
|
@ -150,6 +150,34 @@ private:
|
|||
//!< each incident energy
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Mixed coherent/incoherent elastic angle-energy distribution
|
||||
//==============================================================================
|
||||
|
||||
class MixedElasticAE : public AngleEnergy {
|
||||
public:
|
||||
//! Construct from HDF5 file
|
||||
//
|
||||
//! \param[in] group HDF5 group
|
||||
explicit MixedElasticAE(
|
||||
hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs);
|
||||
|
||||
//! Sample distribution for an angle and energy
|
||||
//! \param[in] E_in Incoming energy in [eV]
|
||||
//! \param[out] E_out Outgoing energy in [eV]
|
||||
//! \param[out] mu Outgoing cosine with respect to current direction
|
||||
//! \param[inout] seed Pseudorandom number seed pointer
|
||||
void sample(
|
||||
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
|
||||
|
||||
private:
|
||||
CoherentElasticAE coherent_dist_; //!< Coherent distribution
|
||||
unique_ptr<AngleEnergy> incoherent_dist_; //!< Incoherent distribution
|
||||
|
||||
const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS
|
||||
const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_SECONDARY_THERMAL_H
|
||||
|
|
|
|||
|
|
@ -44,6 +44,8 @@ class AngleEnergy(EqualityMixin, ABC):
|
|||
return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group)
|
||||
elif dist_type == 'incoherent_inelastic':
|
||||
return openmc.data.IncoherentInelasticAE.from_hdf5(group)
|
||||
elif dist_type == 'mixed_elastic':
|
||||
return openmc.data.MixedElasticAE.from_hdf5(group)
|
||||
|
||||
@staticmethod
|
||||
def from_ace(ace, location_dist, location_start, rx=None):
|
||||
|
|
|
|||
|
|
@ -544,7 +544,7 @@ class Combination(EqualityMixin):
|
|||
self._operations = operations
|
||||
|
||||
|
||||
class Sum(EqualityMixin):
|
||||
class Sum(Function1D):
|
||||
"""Sum of multiple functions.
|
||||
|
||||
This class allows you to create a callable object which represents the sum
|
||||
|
|
@ -578,6 +578,49 @@ class Sum(EqualityMixin):
|
|||
cv.check_type('functions', functions, Iterable, Callable)
|
||||
self._functions = functions
|
||||
|
||||
def to_hdf5(self, group, name='xy'):
|
||||
"""Write sum of functions to an HDF5 group
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to write to
|
||||
name : str
|
||||
Name of the dataset to create
|
||||
|
||||
"""
|
||||
sum_group = group.create_group(name)
|
||||
sum_group.attrs['type'] = np.string_(type(self).__name__)
|
||||
sum_group.attrs['n'] = len(self.functions)
|
||||
for i, f in enumerate(self.functions):
|
||||
f.to_hdf5(sum_group, f'func_{i+1}')
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate sum of functions from an HDF5 group
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
Group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.Sum
|
||||
Functions read from the group
|
||||
|
||||
"""
|
||||
n = group.attrs['n']
|
||||
functions = [
|
||||
Function1D.from_hdf5(group[f'func_{i+1}'])
|
||||
for i in range(n)
|
||||
]
|
||||
return cls(functions)
|
||||
|
||||
|
||||
class Regions1D(EqualityMixin):
|
||||
r"""Piecewise composition of multiple functions.
|
||||
|
|
|
|||
|
|
@ -19,12 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf
|
|||
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, isotopes
|
||||
from .ace import Table, get_table, Library
|
||||
from .angle_energy import AngleEnergy
|
||||
from .function import Tabulated1D, Function1D
|
||||
from .function import Tabulated1D, Function1D, Sum
|
||||
from .njoy import make_ace_thermal
|
||||
from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE,
|
||||
IncoherentElasticAEDiscrete,
|
||||
IncoherentInelasticAEDiscrete,
|
||||
IncoherentInelasticAE)
|
||||
IncoherentInelasticAE, MixedElasticAE)
|
||||
|
||||
|
||||
_THERMAL_NAMES = {
|
||||
|
|
@ -694,29 +694,53 @@ class ThermalScattering(EqualityMixin):
|
|||
|
||||
# Incoherent/coherent elastic scattering cross section
|
||||
idx = ace.jxs[4]
|
||||
n_mu = ace.nxs[6] + 1
|
||||
if idx != 0:
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
|
||||
if ace.nxs[5] == 4:
|
||||
if ace.nxs[5] in (4, 5):
|
||||
# Coherent elastic
|
||||
xs = CoherentElastic(energy, P*EV_PER_MEV)
|
||||
distribution = CoherentElasticAE(xs)
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
coherent_xs = CoherentElastic(energy, P*EV_PER_MEV)
|
||||
coherent_dist = CoherentElasticAE(coherent_xs)
|
||||
|
||||
# Coherent elastic shouldn't have angular distributions listed
|
||||
n_mu = ace.nxs[6] + 1
|
||||
assert n_mu == 0
|
||||
else:
|
||||
# Incoherent elastic
|
||||
xs = Tabulated1D(energy, P)
|
||||
|
||||
if ace.nxs[5] in (3, 5):
|
||||
# Incoherent elastic scattering -- first determine if both
|
||||
# incoherent and coherent are present (mixed)
|
||||
mixed = (ace.nxs[5] == 5)
|
||||
|
||||
# Get cross section values
|
||||
idx = ace.jxs[7] if mixed else ace.jxs[4]
|
||||
n_energy = int(ace.xss[idx])
|
||||
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
|
||||
values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
|
||||
|
||||
incoherent_xs = Tabulated1D(energy, values)
|
||||
|
||||
# Angular distribution
|
||||
n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1
|
||||
assert n_mu > 0
|
||||
idx = ace.jxs[6]
|
||||
idx = ace.jxs[9] if mixed else ace.jxs[6]
|
||||
mu_out = ace.xss[idx:idx + n_energy * n_mu]
|
||||
mu_out.shape = (n_energy, n_mu)
|
||||
distribution = IncoherentElasticAEDiscrete(mu_out)
|
||||
incoherent_dist = IncoherentElasticAEDiscrete(mu_out)
|
||||
|
||||
if ace.nxs[5] == 3:
|
||||
xs = incoherent_xs
|
||||
dist = incoherent_dist
|
||||
elif ace.nxs[5] == 4:
|
||||
xs = coherent_xs
|
||||
dist = coherent_dist
|
||||
else:
|
||||
# Create mixed cross section -- note that coherent must come
|
||||
# first due to assumption on C++ side
|
||||
xs = Sum([coherent_xs, incoherent_xs])
|
||||
|
||||
# Create mixed distribution
|
||||
distribution = MixedElasticAE(coherent_dist, incoherent_dist)
|
||||
|
||||
table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution})
|
||||
|
||||
|
|
@ -802,7 +826,7 @@ class ThermalScattering(EqualityMixin):
|
|||
# Replace ACE data with ENDF data
|
||||
rx, rx_endf = data.elastic, data_endf.elastic
|
||||
for t in temperatures:
|
||||
if isinstance(rx_endf.xs[t], IncoherentElastic):
|
||||
if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)):
|
||||
rx.xs[t] = rx_endf.xs[t]
|
||||
rx.distribution[t] = rx_endf.distribution[t]
|
||||
|
||||
|
|
@ -832,20 +856,14 @@ class ThermalScattering(EqualityMixin):
|
|||
# Read coherent/incoherent elastic data
|
||||
elastic = None
|
||||
if (7, 2) in ev.section:
|
||||
xs = {}
|
||||
distribution = {}
|
||||
|
||||
file_obj = StringIO(ev.section[7, 2])
|
||||
lhtr = endf.get_head_record(file_obj)[2]
|
||||
if lhtr == 1:
|
||||
# coherent elastic
|
||||
|
||||
# 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)
|
||||
xs = {strT: CoherentElastic(bragg_edges, S.y)}
|
||||
distribution = {strT: CoherentElasticAE(xs[strT])}
|
||||
|
||||
# Get structure factor for subsequent temperatures
|
||||
|
|
@ -854,15 +872,34 @@ class ThermalScattering(EqualityMixin):
|
|||
strT = _temperature_str(params[0])
|
||||
xs[strT] = CoherentElastic(bragg_edges, S)
|
||||
distribution[strT] = CoherentElasticAE(xs[strT])
|
||||
return xs, distribution
|
||||
|
||||
elif lhtr == 2:
|
||||
# incoherent elastic
|
||||
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)
|
||||
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@ import numpy as np
|
|||
|
||||
from .angle_energy import AngleEnergy
|
||||
from .correlated import CorrelatedAngleEnergy
|
||||
import openmc.data
|
||||
|
||||
|
||||
class CoherentElasticAE(AngleEnergy):
|
||||
|
|
@ -43,7 +44,27 @@ class CoherentElasticAE(AngleEnergy):
|
|||
|
||||
"""
|
||||
group.attrs['type'] = np.string_('coherent_elastic')
|
||||
group['coherent_xs'] = group.parent['xs']
|
||||
self.coherent_xs.to_hdf5(group, 'coherent_xs')
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate coherent elastic distribution from HDF5 data
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.CoherentElasticAE
|
||||
Coherent elastic distribution
|
||||
|
||||
"""
|
||||
coherent_xs = openmc.data.CoherentElastic.from_hdf5(group['coherent_xs'])
|
||||
return cls(coherent_xs)
|
||||
|
||||
|
||||
class IncoherentElasticAE(AngleEnergy):
|
||||
|
|
@ -101,7 +122,7 @@ class IncoherentElasticAE(AngleEnergy):
|
|||
Incoherent elastic distribution
|
||||
|
||||
"""
|
||||
return cls(group['debye_waller'])
|
||||
return cls(group['debye_waller'][()])
|
||||
|
||||
|
||||
class IncoherentElasticAEDiscrete(AngleEnergy):
|
||||
|
|
@ -210,3 +231,62 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
|
|||
|
||||
class IncoherentInelasticAE(CorrelatedAngleEnergy):
|
||||
_name = 'incoherent_inelastic'
|
||||
|
||||
|
||||
class MixedElasticAE(AngleEnergy):
|
||||
"""Secondary distribution for mixed coherent/incoherent thermal elastic
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
coherent : AngleEnergy
|
||||
Secondary distribution for coherent elastic scattering
|
||||
incoherent : AngleEnergy
|
||||
Secondary distribution for incoherent elastic scattering
|
||||
|
||||
Attributes
|
||||
----------
|
||||
coherent : AngleEnergy
|
||||
Secondary distribution for coherent elastic scattering
|
||||
incoherent : AngleEnergy
|
||||
Secondary distribution for incoherent elastic scattering
|
||||
|
||||
"""
|
||||
def __init__(self, coherent, incoherent):
|
||||
self.coherent = coherent
|
||||
self.incoherent = incoherent
|
||||
|
||||
def to_hdf5(self, group):
|
||||
"""Write mixed elastic distribution to an HDF5 group
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to write to
|
||||
|
||||
"""
|
||||
group.attrs['type'] = np.string_('mixed_elastic')
|
||||
coherent_group = group.create_group('coherent')
|
||||
self.coherent.to_hdf5(coherent_group)
|
||||
incoherent_group = group.create_group('incoherent')
|
||||
self.incoherent.to_hdf5(incoherent_group)
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Generate mixed thermal elastic distribution from HDF5 data
|
||||
|
||||
Parameters
|
||||
----------
|
||||
group : h5py.Group
|
||||
HDF5 group to read from
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.data.MixedElasticAE
|
||||
Mixed thermal elastic distribution
|
||||
|
||||
"""
|
||||
coherent = AngleEnergy.from_hdf5(group['coherent'])
|
||||
incoherent = AngleEnergy.from_hdf5(group['incoherent'])
|
||||
return cls(coherent, incoherent)
|
||||
|
|
|
|||
|
|
@ -14,8 +14,11 @@ class EqualityMixin:
|
|||
def __eq__(self, other):
|
||||
if isinstance(other, type(self)):
|
||||
for key, value in self.__dict__.items():
|
||||
if not np.array_equal(value, other.__dict__.get(key)):
|
||||
return False
|
||||
if isinstance(value, np.ndarray):
|
||||
if not np.array_equal(value, other.__dict__.get(key)):
|
||||
return False
|
||||
else:
|
||||
return value == other.__dict__.get(key)
|
||||
else:
|
||||
return False
|
||||
|
||||
|
|
|
|||
44
src/endf.cpp
44
src/endf.cpp
|
|
@ -90,23 +90,25 @@ bool is_inelastic_scatter(int mt)
|
|||
|
||||
unique_ptr<Function1D> read_function(hid_t group, const char* name)
|
||||
{
|
||||
hid_t dset = open_dataset(group, name);
|
||||
hid_t obj_id = open_object(group, name);
|
||||
std::string func_type;
|
||||
read_attribute(dset, "type", func_type);
|
||||
read_attribute(obj_id, "type", func_type);
|
||||
unique_ptr<Function1D> func;
|
||||
if (func_type == "Tabulated1D") {
|
||||
func = make_unique<Tabulated1D>(dset);
|
||||
func = make_unique<Tabulated1D>(obj_id);
|
||||
} else if (func_type == "Polynomial") {
|
||||
func = make_unique<Polynomial>(dset);
|
||||
func = make_unique<Polynomial>(obj_id);
|
||||
} else if (func_type == "CoherentElastic") {
|
||||
func = make_unique<CoherentElasticXS>(dset);
|
||||
func = make_unique<CoherentElasticXS>(obj_id);
|
||||
} else if (func_type == "IncoherentElastic") {
|
||||
func = make_unique<IncoherentElasticXS>(dset);
|
||||
func = make_unique<IncoherentElasticXS>(obj_id);
|
||||
} else if (func_type == "Sum") {
|
||||
func = make_unique<Sum1D>(obj_id);
|
||||
} else {
|
||||
throw std::runtime_error {"Unknown function type " + func_type +
|
||||
" for dataset " + object_name(dset)};
|
||||
" for dataset " + object_name(obj_id)};
|
||||
}
|
||||
close_dataset(dset);
|
||||
close_object(obj_id);
|
||||
return func;
|
||||
}
|
||||
|
||||
|
|
@ -271,4 +273,30 @@ double IncoherentElasticXS::operator()(double E) const
|
|||
return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Sum1D implementation
|
||||
//==============================================================================
|
||||
|
||||
Sum1D::Sum1D(hid_t group)
|
||||
{
|
||||
// Get number of functions
|
||||
int n;
|
||||
read_attribute(group, "n", n);
|
||||
|
||||
// Get each function
|
||||
for (int i = 0; i < n; ++i) {
|
||||
auto dset_name = fmt::format("func_{}", i + 1);
|
||||
functions_.push_back(read_function(group, dset_name.c_str()));
|
||||
}
|
||||
}
|
||||
|
||||
double Sum1D::operator()(double x) const
|
||||
{
|
||||
double result = 0.0;
|
||||
for (auto& func : functions_) {
|
||||
result += (*func)(x);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -118,6 +118,12 @@ void close_group(hid_t group_id)
|
|||
fatal_error("Failed to close group");
|
||||
}
|
||||
|
||||
void close_object(hid_t obj_id)
|
||||
{
|
||||
if (H5Oclose(obj_id) < 0)
|
||||
fatal_error("Failed to close object");
|
||||
}
|
||||
|
||||
int dataset_ndims(hid_t dset)
|
||||
{
|
||||
hid_t dspace = H5Dget_space(dset);
|
||||
|
|
@ -394,6 +400,12 @@ hid_t open_group(hid_t group_id, const char* name)
|
|||
return H5Gopen(group_id, name, H5P_DEFAULT);
|
||||
}
|
||||
|
||||
hid_t open_object(hid_t group_id, const std::string& name)
|
||||
{
|
||||
ensure_exists(group_id, name.c_str());
|
||||
return H5Oopen(group_id, name.c_str(), H5P_DEFAULT);
|
||||
}
|
||||
|
||||
void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer)
|
||||
{
|
||||
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
|
||||
|
|
|
|||
|
|
@ -332,4 +332,40 @@ void IncoherentInelasticAE::sample(
|
|||
mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// MixedElasticAE implementation
|
||||
//==============================================================================
|
||||
|
||||
MixedElasticAE::MixedElasticAE(
|
||||
hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)
|
||||
: coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)
|
||||
{
|
||||
// Read incoherent elastic distribution
|
||||
hid_t incoherent_group = open_group(group, "incoherent");
|
||||
std::string temp;
|
||||
read_attribute(incoherent_group, "type", temp);
|
||||
if (temp == "incoherent_elastic") {
|
||||
incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);
|
||||
incoherent_dist_ =
|
||||
make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());
|
||||
}
|
||||
close_group(incoherent_group);
|
||||
}
|
||||
|
||||
void MixedElasticAE::sample(
|
||||
double E_in, double& E_out, double& mu, uint64_t* seed) const
|
||||
{
|
||||
// Evaluate coherent and incoherent elastic cross sections
|
||||
double xs_coh = coherent_xs_(E_in);
|
||||
double xs_incoh = incoherent_xs_(E_in);
|
||||
|
||||
if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
|
||||
coherent_dist_.sample(E_in, E_out, mu, seed);
|
||||
} else {
|
||||
incoherent_dist_->sample(E_in, E_out, mu, seed);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -210,14 +210,22 @@ ThermalData::ThermalData(hid_t group)
|
|||
if (temp == "coherent_elastic") {
|
||||
auto xs = dynamic_cast<CoherentElasticXS*>(elastic_.xs.get());
|
||||
elastic_.distribution = make_unique<CoherentElasticAE>(*xs);
|
||||
} else {
|
||||
if (temp == "incoherent_elastic") {
|
||||
elastic_.distribution = make_unique<IncoherentElasticAE>(dgroup);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
|
||||
elastic_.distribution =
|
||||
make_unique<IncoherentElasticAEDiscrete>(dgroup, xs->x());
|
||||
}
|
||||
} else if (temp == "incoherent_elastic") {
|
||||
elastic_.distribution = make_unique<IncoherentElasticAE>(dgroup);
|
||||
} else if (temp == "incoherent_elastic_discrete") {
|
||||
auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
|
||||
elastic_.distribution =
|
||||
make_unique<IncoherentElasticAEDiscrete>(dgroup, xs->x());
|
||||
} else if (temp == "mixed_elastic") {
|
||||
// Get coherent/incoherent cross sections
|
||||
auto mixed_xs = dynamic_cast<Sum1D*>(elastic_.xs.get());
|
||||
const auto& coh_xs =
|
||||
dynamic_cast<const CoherentElasticXS*>(mixed_xs->functions(0).get());
|
||||
const auto& incoh_xs = mixed_xs->functions(1).get();
|
||||
|
||||
// Create mixed elastic distribution
|
||||
elastic_.distribution =
|
||||
make_unique<MixedElasticAE>(dgroup, *coh_xs, *incoh_xs);
|
||||
}
|
||||
|
||||
close_group(elastic_group);
|
||||
|
|
|
|||
|
|
@ -262,3 +262,108 @@ def test_get_thermal_name():
|
|||
|
||||
# Names that don't remotely match anything
|
||||
assert f('boogie_monster') == 'c_boogie_monster'
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def fake_mixed_elastic():
|
||||
fake_tsl = openmc.data.ThermalScattering("c_D_in_7LiD", 1.9968, 4.9, [0.0253])
|
||||
fake_tsl.nuclides = ['H2']
|
||||
|
||||
# Create elastic reaction
|
||||
bragg_edges = [0.00370672, 0.00494229, 0.00988458, 0.01359131, 0.01482688,
|
||||
0.01976918, 0.02347589, 0.02471147, 0.02965376, 0.03336048,
|
||||
0.03953834, 0.04324506, 0.04448063, 0.04942292, 0.05312964,
|
||||
0.05436522, 0.05930751, 0.06301423, 0.0642498 , 0.06919209,
|
||||
0.07289881, 0.07907667, 0.08278339, 0.08401896, 0.08896126,
|
||||
0.09266798, 0.09390355, 0.09884584, 0.1025526 , 0.1037882 ,
|
||||
0.1087305 , 0.1124372 , 0.1186151 , 0.1223218 , 0.1235574 ,
|
||||
0.1284997 , 0.1322064 , 0.133442 , 0.142091 , 0.1433266 ,
|
||||
0.1482688 , 0.1519756 , 0.1581534 , 0.1618601 , 0.1630957 ,
|
||||
0.168038 , 0.1717447 , 0.1729803 , 0.1779226 , 0.1816293 ,
|
||||
0.1828649 , 0.1878072 , 0.1915139 , 0.1976918 , 0.2026341 ,
|
||||
0.2075763 , 0.2125186 , 0.2174609 , 0.2224032 , 0.2273455 ,
|
||||
0.2421724 , 0.2471147 , 0.252057 , 0.2569993 , 0.2619415 ,
|
||||
0.2668838 , 0.2767684 , 0.2817107 , 0.2915953 , 0.3064222 ,
|
||||
0.3261913 , 0.366965]
|
||||
factors = [0.00375735, 0.01386287, 0.02595574, 0.02992438, 0.03549502,
|
||||
0.03855745, 0.04058831, 0.04986305, 0.05703106, 0.05855471,
|
||||
0.06078031, 0.06212291, 0.06656602, 0.06930339, 0.0697072 ,
|
||||
0.07201456, 0.07263853, 0.07313129, 0.07465531, 0.07714482,
|
||||
0.07759976, 0.077809 , 0.07790282, 0.07927957, 0.08013058,
|
||||
0.08026637, 0.08073475, 0.08112202, 0.08123039, 0.08187171,
|
||||
0.08213756, 0.08218236, 0.08236572, 0.08240729, 0.08259795,
|
||||
0.08297893, 0.08300455, 0.08314566, 0.08315611, 0.08337715,
|
||||
0.08350026, 0.08350663, 0.08352815, 0.08353776, 0.0836098 ,
|
||||
0.08367017, 0.08367361, 0.0837242 , 0.08375069, 0.08375227,
|
||||
0.08377006, 0.08381488, 0.08381644, 0.08382698, 0.08386266,
|
||||
0.08387756, 0.08388445, 0.08388974, 0.08390341, 0.08391088,
|
||||
0.08391695, 0.08392361, 0.08392684, 0.08392818, 0.08393161,
|
||||
0.08393546, 0.08393685, 0.08393801, 0.08393976, 0.08394167,
|
||||
0.08394288, 0.08394398]
|
||||
coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors)
|
||||
incoherent_xs = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672])
|
||||
elastic_xs = {'294K': openmc.data.Sum((coherent_xs, incoherent_xs))}
|
||||
coherent_dist = openmc.data.CoherentElasticAE(coherent_xs)
|
||||
incoherent_dist = openmc.data.IncoherentElasticAEDiscrete([
|
||||
[-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6]
|
||||
])
|
||||
elastic_dist = {'294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist)}
|
||||
fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist)
|
||||
|
||||
# Create inelastic reaction
|
||||
inelastic_xs = {'294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35])}
|
||||
breakpoints = [3]
|
||||
interpolation = [2]
|
||||
energy = [1.0e-5, 4.3e-2, 4.9]
|
||||
energy_out = [
|
||||
openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]),
|
||||
openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]),
|
||||
openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]),
|
||||
]
|
||||
for eout in energy_out:
|
||||
eout.normalize()
|
||||
eout.c = eout.cdf()
|
||||
discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8)
|
||||
discrete.c = discrete.cdf()[1:]
|
||||
mu = [[discrete]*4]*3
|
||||
inelastic_dist = {'294K': openmc.data.IncoherentInelasticAE(
|
||||
breakpoints, interpolation, energy, energy_out, mu)}
|
||||
inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist)
|
||||
fake_tsl.inelastic = inelastic
|
||||
|
||||
return fake_tsl
|
||||
|
||||
|
||||
def test_mixed_elastic(fake_mixed_elastic, run_in_tmpdir):
|
||||
# Write data to HDF5 and then read back
|
||||
original = fake_mixed_elastic
|
||||
original.export_to_hdf5('c_D_in_7LiD.h5')
|
||||
copy = openmc.data.ThermalScattering.from_hdf5('c_D_in_7LiD.h5')
|
||||
|
||||
# Make sure data did not change as a result of HDF5 writing/reading
|
||||
assert original == copy
|
||||
|
||||
# Create modified cross_sections.xml file that includes the above data
|
||||
xs = openmc.data.DataLibrary.from_xml()
|
||||
xs.register_file('c_D_in_7LiD.h5')
|
||||
xs.export_to_xml('cross_sections_mixed.xml')
|
||||
|
||||
# Create a minimal model that includes the new data and run it
|
||||
mat = openmc.Material()
|
||||
mat.add_nuclide('H2', 1.0)
|
||||
mat.add_nuclide('Li7', 1.0)
|
||||
mat.set_density('g/cm3', 1.0)
|
||||
mat.add_s_alpha_beta('c_D_in_7LiD')
|
||||
sph = openmc.Sphere(r=10.0, boundary_type="vacuum")
|
||||
cell = openmc.Cell(fill=mat, region=-sph)
|
||||
model = openmc.Model()
|
||||
model.geometry = openmc.Geometry([cell])
|
||||
model.materials = openmc.Materials([mat])
|
||||
model.materials.cross_sections = "cross_sections_mixed.xml"
|
||||
model.settings.particles = 1000
|
||||
model.settings.batches = 10
|
||||
model.settings.run_mode = 'fixed source'
|
||||
model.settings.source = openmc.Source(
|
||||
energy=openmc.stats.Discrete([3.0], [1.0]) # 3 eV source
|
||||
)
|
||||
model.run()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue