crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into dep-timings

Browse source
This commit is contained in:
Andrew Johnson 2019-07-02 14:02:18 -05:00
commit 975695b9bc
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
36 changed files with 1639 additions and 628 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -219,6 +219,7 @@ list(APPEND libopenmc_SOURCES
src/secondary_correlated.cpp
src/secondary_kalbach.cpp
src/secondary_nbody.cpp
src/secondary_thermal.cpp
src/secondary_uncorrelated.cpp
src/settings.cpp
src/simulation.cpp

107
docs/source/io_formats/nuclear_data.rst
View file

@ -215,11 +215,9 @@ Thermal Neutron Scattering Data
**/<thermal name>/**
:Attributes: - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses
- **nuclides** (*char[][]*) -- Names of nuclides for which the thermal
scattering data applies to
- **secondary_mode** (*char[]*) -- Indicates how the inelastic
outgoing angle-energy distributions are represented ('equal',
'skewed', or 'continuous').
- **energy_max** (*double*) -- Maximum energy in [eV]
- **nuclides** (*char[][]*) -- Names of nuclides for which the
thermal scattering data applies to
**/<thermal name>/kTs/**
@ -237,11 +235,13 @@ temperature-dependent data set. For example, the data set corresponding to
temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Datasets: - **xs** (:ref:`tabulated <1d_tabulated>`) -- Thermal inelastic
:Datasets:
- **xs** (:ref:`function <1d_functions>`) -- Thermal elastic
scattering cross section for temperature TTT (in Kelvin)
- **mu_out** (*double[][]*) -- Distribution of outgoing energies
and angles for coherent elastic scattering for temperature TTT
(in Kelvin)
:Groups:
- **distribution** -- Format for angle-energy distributions are
detailed in :ref:`angle_energy`.
**/<thermal name>/inelastic/<TTT>K/**
@ -249,18 +249,13 @@ temperature-dependent data set. For example, the data set corresponding to
temperature-dependent data set. For example, the data set corresponding to
300 Kelvin would be located at `300K`.
:Datasets: - **xs** (:ref:`tabulated <1d_tabulated>`) -- Thermal inelastic
:Datasets:
- **xs** (:ref:`function <1d_functions>`) -- Thermal inelastic
scattering cross section for temperature TTT (in Kelvin)
- **energy_out** (*double[][]*) -- Distribution of outgoing
energies for each incoming energy for temperature TTT (in Kelvin).
Only present if secondary mode is not continuous.
- **mu_out** (*double[][][]*) -- Distribution of scattering cosines
for each pair of incoming and outgoing energies. for temperature
TTT (in Kelvin). Only present if secondary mode is not continuous.
If the secondary mode is continuous, the outgoing energy-angle distribution is
given as a :ref:`correlated angle-energy distribution
<correlated_angle_energy>`.
:Groups:
- **distribution** -- Format for angle-energy distributions are
detailed in :ref:`angle_energy`.
.. _product:
@ -332,7 +327,17 @@ Coherent elastic scattering
:Datatype: *double[2][]*
:Description: The first row lists Bragg edges and the second row lists structure
factor cumulative sums.
:Attributes: - **type** (*char[]*) -- 'bragg'
:Attributes: - **type** (*char[]*) -- 'CoherentElastic'
Incoherent elastic scattering
-----------------------------
:Object type: Dataset
:Datatype: *double[2]*
:Description: The first value is the characteristic bound cross section in [b]
and the second value is the Debye-Waller integral in
[eV\ :math:`^{-1}`].
:Attributes: - **type** (*char[]*) -- 'IncoherentElastic'
.. _angle_energy:
@ -434,6 +439,68 @@ N-Body Phase Space
target nuclide in neutron masses
- **q_value** (*double*) -- Q value for the reaction in eV
Coherent Elastic
----------------
This angle-energy distribution is used specifically for coherent elastic thermal
neutron scattering.
:Object type: Group
:Attributes: - **type** (*char[]*) -- "coherent_elastic"
:Hard link: - **xs** -- Link to the coherent elastic scattering cross section
Incoherent Elastic
------------------
This angle-energy distribution is used specifically for incoherent elastic
thermal neutron scattering (derived from an ENDF file directly).
:Object type: Group
:Attributes: - **type** (*char[]*) -- "incoherent_elastic"
:Datasets:
- **debye_waller** (*double*) -- Debye-Waller integral in
[eV\ :math:`^{-1}`]
Incoherent Elastic (Discrete)
-----------------------------
This angle-energy distribution is used for discretized incoherent elastic
thermal neutron scattering distributions that are present in ACE files.
:Object type: Group
:Attributes: - **type** (*char[]*) -- "incoherent_elastic_discrete"
:Datasets:
- **mu_out** (*double[][]*) -- Equiprobable discrete outgoing
angles for each incident neutron energy tabulated
Incoherent Inelastic
--------------------
This angle-energy distribution is used specifically for (continuous) incoherent
inelastic thermal neutron scattering.
:Object type: Group
:Attributes: - **type** (*char[]*) -- "incoherent_inelastic"
:Datasets: The datasets for this angle-energy distribution are the same as for
:ref:`correlated angle-energy distributions
<correlated_angle_energy>`.
Incoherent Inelastic (Discrete)
-------------------------------
This angle-energy distribution is used specifically for incoherent inelastic
thermal neutron scattering where the distributions have been discretized into
equiprobable bins.
:Object type: Group
:Attributes: - **type** (*char[]*) -- "incoherent_inelastic_discrete"
:Datasets: - **energy_out** (*double[][]*) -- Distribution of outgoing
energies for each incoming energy.
- **mu_out** (*double[][][]*) -- Distribution of scattering cosines
for each pair of incoming and outgoing energies.
- **skewed** (*int8_t*) -- Whether discrete angles are equi-probable
(0) or have a skewed distribution (1).
.. _energy_distribution:
--------------------

57
docs/source/pythonapi/data.rst
View file

@ -5,6 +5,9 @@
Core Classes
------------
The following classes are used for incident neutron data, decay data, fission
and product yields.
.. autosummary::
:toctree: generated
:nosignatures:
@ -13,18 +16,38 @@ Core Classes
openmc.data.IncidentNeutron
openmc.data.Reaction
openmc.data.Product
openmc.data.Tabulated1D
openmc.data.FissionEnergyRelease
openmc.data.ThermalScattering
openmc.data.CoherentElastic
openmc.data.FissionEnergyRelease
openmc.data.DataLibrary
openmc.data.IncidentPhoton
openmc.data.PhotonReaction
openmc.data.AtomicRelaxation
openmc.data.Decay
openmc.data.FissionProductYields
openmc.data.WindowedMultipole
openmc.data.ProbabilityTables
The following classes are used for storing atomic data (incident photon cross
sections, atomic relaxation):
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.data.IncidentPhoton
openmc.data.PhotonReaction
openmc.data.AtomicRelaxation
The following classes are used for storing thermal neutron scattering data:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.data.ThermalScattering
openmc.data.ThermalScatteringReaction
openmc.data.CoherentElastic
openmc.data.IncoherentElastic
Core Functions
--------------
@ -41,6 +64,22 @@ Core Functions
openmc.data.water_density
openmc.data.zam
One-dimensional Functions
-------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.data.Function1D
openmc.data.Tabulated1D
openmc.data.Polynomial
openmc.data.Combination
openmc.data.Sum
openmc.data.Regions1D
openmc.data.ResonancesWithBackground
Angle-Energy Distributions
--------------------------
@ -66,6 +105,10 @@ Angle-Energy Distributions
openmc.data.DiscretePhoton
openmc.data.LevelInelastic
openmc.data.ContinuousTabular
openmc.data.CoherentElasticAE
openmc.data.IncoherentElasticAE
openmc.data.IncoherentElasticAEDiscrete
openmc.data.IncoherentInelasticAEDiscrete
Resonance Data
--------------

8
docs/source/usersguide/cross_sections.rst
View file

@ -11,7 +11,7 @@ or multi-group mode.
In continuous-energy mode, OpenMC uses a native `HDF5
<https://support.hdfgroup.org/HDF5/>`_ format (see :ref:`io_nuclear_data`) to
store all nuclear data. Pregenerated HDF5 libraries can be found at
https://openmc.mcs.anl.gov; unless you have specific data needs, it is highly
https://openmc.org; unless you have specific data needs, it is highly
recommended to use one of the pregenerated libraries. Alternatively, if you have
ACE format data that was produced with NJOY_, such as that distributed with
MCNP_ or Serpent_, it can be converted to the HDF5 format using the :ref:`using
@ -61,7 +61,7 @@ Using Pregenerated Libraries
----------------------------
Various evaluated nuclear data libraries have been processed into the HDF5
format required by OpenMC and can be found at https://openmc.mcs.anl.gov. You
format required by OpenMC and can be found at https://openmc.org. You
can find both libraries generated by the OpenMC development team as well as
libraries based on ACE files distributed elsewhere. To use these libraries,
download the archive file, unpack it, and then set your
@ -198,7 +198,7 @@ file, is distributed with OpenMC. The rest is available from the NNDC_, which
provides ENDF data from the photo-atomic and atomic relaxation sublibraries of
the ENDF/B-VII.1 library.
Most of the pregenerated HDF5 libraries available at https://openmc.mcs.anl.gov
Most of the pregenerated HDF5 libraries available at https://openmc.org
already have photon interaction data included. If you are building a data
library yourself, it is possible to use the Python API directly to convert
photon interaction data from an ENDF or ACE file to an HDF5 file. The
@ -237,7 +237,7 @@ the :class:`openmc.data.DataLibrary` class to register the .h5 files as
described in :ref:`create_xs_library`.
The `official ENDF/B-VII.1 HDF5 library
<https://openmc.mcs.anl.gov/official-data-libraries/>`_ includes the windowed
<https://openmc.org/official-data-libraries/>`_ includes the windowed
multipole library, so if you are using this library, the windowed multipole data
will already be available to you.

3
include/openmc/constants.h
View file

@ -26,7 +26,7 @@ constexpr bool VERSION_DEV {true};
constexpr std::array<int, 3> VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE};
// HDF5 data format
constexpr int HDF5_VERSION[] {2, 0};
constexpr int HDF5_VERSION[] {3, 0};
// Version numbers for binary files
constexpr std::array<int, 2> VERSION_STATEPOINT {17, 0};
@ -36,7 +36,6 @@ constexpr std::array<int, 2> VERSION_SUMMARY {6, 0};
constexpr std::array<int, 2> VERSION_VOLUME {1, 0};
constexpr std::array<int, 2> VERSION_VOXEL {2, 0};
constexpr std::array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
constexpr char VERSION_MULTIPOLE[] {"v0.2"};
// ============================================================================
// ADJUSTABLE PARAMETERS

25
include/openmc/endf.h
View file

@ -57,7 +57,7 @@ public:
//! Evaluate the polynomials
//! \param[in] x independent variable
//! \return Polynomial evaluated at x
double operator()(double x) const;
double operator()(double x) const override;
private:
std::vector<double> coef_; //!< Polynomial coefficients
};
@ -77,7 +77,7 @@ public:
//! Evaluate the tabulated function
//! \param[in] x independent variable
//! \return Function evaluated at x
double operator()(double x) const;
double operator()(double x) const override;
// Accessors
const std::vector<double>& x() const { return x_; }
@ -96,13 +96,32 @@ private:
//==============================================================================
class CoherentElasticXS : public Function1D {
public:
explicit CoherentElasticXS(hid_t dset);
double operator()(double E) const;
double operator()(double E) const override;
const std::vector<double>& bragg_edges() const { return bragg_edges_; }
const std::vector<double>& factors() const { return factors_; }
private:
std::vector<double> bragg_edges_; //!< Bragg edges in [eV]
std::vector<double> factors_; //!< Partial sums of structure factors [eV-b]
};
//==============================================================================
//! Incoherent elastic scattering cross section
//==============================================================================
class IncoherentElasticXS : public Function1D {
public:
explicit IncoherentElasticXS(hid_t dset);
double operator()(double E) const override;
private:
double bound_xs_; //!< Characteristic bound xs in [b]
double debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1]
};
//! Read 1D function from HDF5 dataset
//! \param[in] group HDF5 group containing dataset
//! \param[in] name Name of dataset

2
include/openmc/secondary_correlated.h
View file

@ -38,7 +38,7 @@ public:
//! \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
void sample(double E_in, double& E_out, double& mu) const;
void sample(double E_in, double& E_out, double& mu) const override;
// energy property
std::vector<double>& energy() { return energy_; }

2
include/openmc/secondary_kalbach.h
View file

@ -29,7 +29,7 @@ public:
//! \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
void sample(double E_in, double& E_out, double& mu) const;
void sample(double E_in, double& E_out, double& mu) const override;
private:
//! Outgoing energy/angle at a single incoming energy
struct KMTable {

2
include/openmc/secondary_nbody.h
View file

@ -24,7 +24,7 @@ public:
//! \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
void sample(double E_in, double& E_out, double& mu) const;
void sample(double E_in, double& E_out, double& mu) const override;
private:
int n_bodies_; //!< Number of particles distributed
double mass_ratio_; //!< Total mass of particles [neutron mass]

139
include/openmc/secondary_thermal.h Normal file
View file

@ -0,0 +1,139 @@
//! \file secondary_thermal.h
//! Angle-energy distributions for thermal scattering
#ifndef OPENMC_SECONDARY_THERMAL_H
#define OPENMC_SECONDARY_THERMAL_H
#include "openmc/angle_energy.h"
#include "openmc/endf.h"
#include "openmc/secondary_correlated.h"
#include <hdf5.h>
#include "xtensor/xtensor.hpp"
#include <vector>
namespace openmc {
//==============================================================================
//! Coherent elastic scattering angle-energy distribution
//==============================================================================
class CoherentElasticAE : public AngleEnergy {
public:
//! Construct from a coherent elastic scattering cross section
//
//! \param[in] xs Coherent elastic scattering cross section
explicit CoherentElasticAE(const CoherentElasticXS& 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
void sample(double E_in, double& E_out, double& mu) const override;
private:
const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
};
//==============================================================================
//! Incoherent elastic scattering angle-energy distribution
//==============================================================================
class IncoherentElasticAE : public AngleEnergy {
public:
//! Construct from HDF5 file
//
//! \param[in] group HDF5 group
explicit IncoherentElasticAE(hid_t group);
//! 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
void sample(double E_in, double& E_out, double& mu) const override;
private:
double debye_waller_;
};
//==============================================================================
//! Incoherent elastic scattering angle-energy distribution (discrete)
//==============================================================================
class IncoherentElasticAEDiscrete : public AngleEnergy {
public:
//! Construct from HDF5 file
//
//! \param[in] group HDF5 group
//! \param[in] energy Energies at which cosines are tabulated
explicit IncoherentElasticAEDiscrete(hid_t group, const std::vector<double>& energy);
//! 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
void sample(double E_in, double& E_out, double& mu) const override;
private:
const std::vector<double>& energy_; //!< Energies at which cosines are tabulated
xt::xtensor<double, 2> mu_out_; //!< Cosines for each incident energy
};
//==============================================================================
//! Incoherent inelastic scattering angle-energy distribution (discrete)
//==============================================================================
class IncoherentInelasticAEDiscrete : public AngleEnergy {
public:
//! Construct from HDF5 file
//
//! \param[in] group HDF5 group
//! \param[in] energy Incident energies at which distributions are tabulated
explicit IncoherentInelasticAEDiscrete(hid_t group, const std::vector<double>& energy);
//! 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
void sample(double E_in, double& E_out, double& mu) const override;
private:
const std::vector<double>& energy_; //!< Incident energies
xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
xt::xtensor<double, 3> mu_out_; //!< Outgoing cosines for each incident/outgoing energy
bool skewed_; //!< Whether outgoing energy distribution is skewed
};
//==============================================================================
//! Incoherent inelastic scattering angle-energy distribution
//==============================================================================
class IncoherentInelasticAE : public AngleEnergy {
public:
//! Construct from HDF5 file
//
//! \param[in] group HDF5 group
explicit IncoherentInelasticAE(hid_t group);
//! 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
void sample(double E_in, double& E_out, double& mu) const override;
private:
//! Secondary energy/angle distribution
struct DistEnergySab {
std::size_t n_e_out; //!< Number of outgoing energies
xt::xtensor<double, 1> e_out; //!< Outgoing energies
xt::xtensor<double, 1> e_out_pdf; //!< Probability density function
xt::xtensor<double, 1> e_out_cdf; //!< Cumulative distribution function
xt::xtensor<double, 2> mu; //!< Equiprobable angles at each outgoing energy
};
std::vector<double> energy_; //!< Incident energies
std::vector<DistEnergySab> distribution_; //!< Secondary angle-energy at
//!< each incident energy
};
} // namespace openmc
#endif // OPENMC_SECONDARY_THERMAL_H

2
include/openmc/secondary_uncorrelated.h
View file

@ -29,7 +29,7 @@ public:
//! \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
void sample(double E_in, double& E_out, double& mu) const;
void sample(double E_in, double& E_out, double& mu) const override;
// Accessors
AngleDistribution& angle() { return angle_; }

62
include/openmc/thermal.h
View file

@ -9,6 +9,8 @@
#include "xtensor/xtensor.hpp"
#include "openmc/angle_energy.h"
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
@ -47,20 +49,31 @@ extern std::unordered_map<std::string, int> thermal_scatt_map;
class ThermalData {
public:
ThermalData(hid_t group, int secondary_mode);
ThermalData(hid_t group);
// Sample an outgoing energy and angle
//! Calculate the cross section
//
//! \param[in] E Incident neutron energy in [eV]
//! \param[out] elastic Elastic scattering cross section in [b]
//! \param[out] inelastic Inelastic scattering cross section in [b]
void calculate_xs(double E, double* elastic, double* inelastic) const;
//! Sample an outgoing energy and angle
//
//! \param[in] micro_xs Microscopic cross sections
//! \param[in] E_in Incident neutron energy in [eV]
//! \param[out] E_out Outgoing neutron energy in [eV]
//! \param[out] mu Outgoing scattering angle cosine
void sample(const NuclideMicroXS& micro_xs, double E_in,
double* E_out, double* mu);
private:
//! Secondary energy/angle distributions for inelastic thermal scattering
//! collisions which utilize a continuous secondary energy representation.
struct DistEnergySab {
std::size_t n_e_out; //!< Number of outgoing energies
xt::xtensor<double, 1> e_out; //!< Outgoing energies
xt::xtensor<double, 1> e_out_pdf; //!< Probability density function
xt::xtensor<double, 1> e_out_cdf; //!< Cumulative distribution function
xt::xtensor<double, 2> mu; //!< Equiprobable angles at each outgoing energy
struct Reaction {
// Default constructor
Reaction() { }
// Data members
std::unique_ptr<Function1D> xs; //!< Cross section
std::unique_ptr<AngleEnergy> distribution; //!< Secondary angle-energy distribution
};
//! Upper threshold for incoherent inelastic scattering (usually ~4 eV)
@ -69,30 +82,8 @@ private:
double threshold_elastic_ {0.0};
// Inelastic scattering data
int inelastic_mode_; //!< distribution type (equal/skewed/continuous)
std::size_t n_inelastic_e_in_; //!< number of incoming E for inelastic
std::size_t n_inelastic_e_out_; //!< number of outgoing E for inelastic
std::size_t n_inelastic_mu_; //!< number of outgoing angles for inelastic
std::vector<double> inelastic_e_in_; //!< incoming E grid for inelastic
std::vector<double> inelastic_sigma_; //!< inelastic scattering cross section
// The following are used only for equal/skewed distributions
xt::xtensor<double, 2> inelastic_e_out_;
xt::xtensor<double, 3> inelastic_mu_;
// The following is used only for continuous S(a,b) distributions. The
// different implementation is necessary because the continuous representation
// has a variable number of outgoing energy points for each incoming energy
std::vector<DistEnergySab> inelastic_data_; //!< Secondary angle-energy at
//!< each incoming energy
// Elastic scattering data
int elastic_mode_; //!< type of elastic (incoherent/coherent)
std::size_t n_elastic_e_in_; //!< number of incoming E for elastic
std::size_t n_elastic_mu_; //!< number of outgoing angles for elastic
std::vector<double> elastic_e_in_; //!< incoming E grid for elastic
std::vector<double> elastic_P_; //!< elastic scattering cross section
xt::xtensor<double, 2> elastic_mu_; //!< equi-probable angles at each incoming E
Reaction elastic_;
Reaction inelastic_;
// ThermalScattering needs access to private data members
friend class ThermalScattering;
@ -127,10 +118,9 @@ public:
void sample(const NuclideMicroXS& micro_xs, double E_in,
double* E_out, double* mu);
double threshold() const { return data_[0].threshold_inelastic_; }
std::string name_; //!< name of table, e.g. "c_H_in_H2O"
double awr_; //!< weight of nucleus in neutron masses
double energy_max_; //!< maximum energy for thermal scattering in [eV]
std::vector<double> kTs_; //!< temperatures in [eV] (k*T)
std::vector<std::string> nuclides_; //!< Valid nuclides

4
openmc/data/__init__.py
View file

@ -1,5 +1,5 @@
# Version of HDF5 nuclear data format
HDF5_VERSION_MAJOR = 2
HDF5_VERSION_MAJOR = 3
HDF5_VERSION_MINOR = 0
HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
@ -16,7 +16,6 @@ from .decay import *
from .reaction import *
from . import ace
from .angle_distribution import *
from .function import *
from . import endf
from .energy_distribution import *
from .product import *
@ -33,3 +32,4 @@ from .resonance import *
from .resonance_covariance import *
from .multipole import *
from .grid import *
from .function import *

10
openmc/data/angle_energy.py
View file

@ -35,6 +35,16 @@ class AngleEnergy(EqualityMixin, metaclass=ABCMeta):
return openmc.data.KalbachMann.from_hdf5(group)
elif dist_type == 'nbody':
return openmc.data.NBodyPhaseSpace.from_hdf5(group)
elif dist_type == 'coherent_elastic':
return openmc.data.CoherentElasticAE.from_hdf5(group)
elif dist_type == 'incoherent_elastic':
return openmc.data.IncoherentElasticAE.from_hdf5(group)
elif dist_type == 'incoherent_elastic_discrete':
return openmc.data.IncoherentElasticAEDiscrete.from_hdf5(group)
elif dist_type == 'incoherent_inelastic_discrete':
return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group)
elif dist_type == 'incoherent_inelastic':
return openmc.data.IncoherentInelasticAE.from_hdf5(group)
@staticmethod
def from_ace(ace, location_dist, location_start, rx=None):

4
openmc/data/correlated.py
View file

@ -44,6 +44,8 @@ class CorrelatedAngleEnergy(AngleEnergy):
"""
_name = 'correlated'
def __init__(self, breakpoints, interpolation, energy, energy_out, mu):
super().__init__()
self.breakpoints = breakpoints
@ -111,7 +113,7 @@ class CorrelatedAngleEnergy(AngleEnergy):
HDF5 group to write to
"""
group.attrs['type'] = np.string_('correlated')
group.attrs['type'] = np.string_(self._name)
dset = group.create_dataset('energy', data=self.energy)
dset.attrs['interpolation'] = np.vstack((self.breakpoints,

3
openmc/data/decay.py
View file

@ -423,7 +423,8 @@ class Decay(EqualityMixin):
if spectrum['type'] == 'ec/beta+':
di['positron_intensity'] = ufloat(*values[4:6])
elif spectrum['type'] == 'gamma':
di['internal_pair'] = ufloat(*values[4:6])
if len(values) >= 6:
di['internal_pair'] = ufloat(*values[4:6])
if len(values) >= 8:
di['total_internal_conversion'] = ufloat(*values[6:8])
if len(values) == 12:

12
openmc/data/function.py
View file

@ -440,6 +440,14 @@ class Tabulated1D(Function1D):
class Polynomial(np.polynomial.Polynomial, Function1D):
"""A power series class.
Parameters
----------
coef : Iterable of float
Polynomial coefficients in order of increasing degree
"""
def to_hdf5(self, group, name='xy'):
"""Write polynomial function to an HDF5 group
@ -583,8 +591,8 @@ class Regions1D(EqualityMixin):
Functions which are to be combined in a piecewise fashion
breakpoints : Iterable of float
The values of the dependent variable that define the domain of
each function. The *i*th and *(i+1)*th values are the limits of the
domain of the *i*th function. Values must be monotonically increasing.
each function. The `i`\ th and `(i+1)`\ th values are the limits of the
domain of the `i`\ th function. Values must be monotonically increasing.
Attributes
----------

8
openmc/data/njoy.py
View file

@ -120,7 +120,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
'{library}: {zsymam_thermal} processed by NJOY'/
{mat} {temperature} '{data.name}' /
{zaids} /
222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} 2/
222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/
"""
@ -355,8 +355,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None,
def make_ace_thermal(filename, filename_thermal, temperatures=None,
ace='ace', xsdir='xsdir', error=0.001, evaluation=None,
evaluation_thermal=None, **kwargs):
ace='ace', xsdir='xsdir', error=0.001, iwt=2,
evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate thermal scattering ACE file from ENDF files
Parameters
@ -374,6 +374,8 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
Path of xsdir file to write
error : float, optional
Fractional error tolerance for NJOY processing
iwt : int
`iwt` parameter used in NJOR/ACER card 9
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.

7
openmc/data/reaction.py
View file

@ -894,10 +894,7 @@ class Reaction(EqualityMixin):
Tgroup = group.create_group(T)
if self.xs[T] is not None:
dset = Tgroup.create_dataset('xs', data=self.xs[T].y)
if hasattr(self.xs[T], '_threshold_idx'):
threshold_idx = self.xs[T]._threshold_idx + 1
else:
threshold_idx = 1
threshold_idx = getattr(self.xs[T], '_threshold_idx', 0)
dset.attrs['threshold_idx'] = threshold_idx
for i, p in enumerate(self.products):
pgroup = group.create_group('product_{}'.format(i))
@ -939,7 +936,7 @@ class Reaction(EqualityMixin):
'at T={} because no corresponding energy grid '
'exists.'.format(mt, T))
xs = Tgroup['xs'][()]
threshold_idx = Tgroup['xs'].attrs['threshold_idx'] - 1
threshold_idx = Tgroup['xs'].attrs['threshold_idx']
tabulated_xs = Tabulated1D(energy[T][threshold_idx:], xs)
tabulated_xs._threshold_idx = threshold_idx
rx.xs[T] = tabulated_xs

551
openmc/data/thermal.py
View file

@ -1,6 +1,8 @@
from collections.abc import Iterable
from collections import namedtuple
from difflib import get_close_matches
from numbers import Real
from io import StringIO
import itertools
import os
import re
@ -13,13 +15,16 @@ import h5py
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE
from .ace import Table, get_table, Library
from .angle_energy import AngleEnergy
from .function import Tabulated1D
from .correlated import CorrelatedAngleEnergy
from .function import Tabulated1D, Function1D
from .njoy import make_ace_thermal
from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE,
IncoherentElasticAEDiscrete,
IncoherentInelasticAEDiscrete,
IncoherentInelasticAE)
_THERMAL_NAMES = {
@ -69,6 +74,12 @@ _THERMAL_NAMES = {
}
def _temperature_str(T):
# round() normally returns an int when called with a single argument, but
# numpy floats overload rounding to return another float
return "{}K".format(int(round(T)))
def get_thermal_name(name):
"""Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'
@ -110,31 +121,40 @@ def get_thermal_name(name):
break
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(name, match))
'Assigning a name of {}.'.format(name, match))
return match
else:
# OK, we give up. Just use the ACE name.
warn('Thermal scattering material "{0}" is not recognized. '
'Assigning a name of c_{0}.'.format(name))
'Assigning a name of c_{0}.'.format(name))
return 'c_' + name
class CoherentElastic(EqualityMixin):
class CoherentElastic(Function1D):
r"""Coherent elastic scattering data from a crystalline material
The integrated cross section for coherent elastic scattering from a
powdered crystalline material may be represented as:
.. math::
\sigma(E,T) = \frac{1}{E} \sum\limits_{i=1}^{E_i < E} s_i(T)
where :math:`s_i(T)` is proportional the structure factor in [eV-b] at
the moderator temperature :math:`T` in Kelvin.
Parameters
----------
bragg_edges : Iterable of float
Bragg edge energies in eV
factors : Iterable of float
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} S_i`
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} s_i`
Attributes
----------
bragg_edges : Iterable of float
Bragg edge energies in eV
factors : Iterable of float
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} S_i`
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} s_i`
"""
@ -143,10 +163,15 @@ class CoherentElastic(EqualityMixin):
self.factors = factors
def __call__(self, E):
idx = np.searchsorted(self.bragg_edges, E) - 1
if isinstance(E, Iterable):
E = np.asarray(E)
idx = np.searchsorted(self.bragg_edges, E)
return self.factors[idx] / E
nonzero = idx >= 0
xs = np.zeros_like(E)
xs[nonzero] = self.factors[idx[nonzero]] / E[nonzero]
return xs
else:
return self.factors[idx] / E if idx >= 0 else 0.0
def __len__(self):
return len(self.bragg_edges)
@ -183,7 +208,7 @@ class CoherentElastic(EqualityMixin):
"""
dataset = group.create_dataset(name, data=np.vstack(
[self.bragg_edges, self.factors]))
dataset.attrs['type'] = np.string_('bragg')
dataset.attrs['type'] = np.string_(type(self).__name__)
@classmethod
def from_hdf5(cls, dataset):
@ -191,8 +216,8 @@ class CoherentElastic(EqualityMixin):
Parameters
----------
group : h5py.Dataset
HDF5 group to write to
dataset : h5py.Dataset
HDF5 dataset to read from
Returns
-------
@ -205,6 +230,152 @@ class CoherentElastic(EqualityMixin):
return cls(bragg_edges, factors)
class IncoherentElastic(Function1D):
r"""Incoherent elastic scattering cross section
Elastic scattering can be treated in the incoherent approximation for
partially ordered systems such as ZrHx and polyethylene. The integrated
cross section can be obtained as:
.. math::
\sigma(E,T) = \frac{\sigma_b}{2} \left ( \frac{1 - e^{-4EW'(T)}}
{2EW'(T)} \right )
where :math:`\sigma_b` is the characteristic bound cross section, and
:math:`W'(T)` is the Debye-Waller integral divided by the atomic mass
in [eV\ :math:`^{-1}`].
Parameters
----------
bound_xs : float
Characteristic bound cross section in [b]
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
Attributes
----------
bound_xs : float
Characteristic bound cross section in [b]
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
"""
def __init__(self, bound_xs, debye_waller):
self.bound_xs = bound_xs
self.debye_waller = debye_waller
def __call__(self, E):
W = self.debye_waller
return self.bound_xs / 2.0 * (1 - np.exp(-4*E*W)) / (2*E*W)
def to_hdf5(self, group, name):
"""Write incoherent elastic scattering to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : str
Name of the dataset to create
"""
data = np.array([self.bound_xs, self.debye_waller])
dataset = group.create_dataset(name, data=data)
dataset.attrs['type'] = np.string_(type(self).__name__)
@classmethod
def from_hdf5(cls, dataset):
"""Read incoherent elastic scattering from an HDF5 dataset
Parameters
----------
dataset : h5py.Dataset
HDF5 dataset to read from
Returns
-------
openmc.data.IncoherentElastic
Incoherent elastic scattering cross section
"""
bound_xs, debye_waller = dataset[()]
return cls(bound_xs, debye_waller)
class ThermalScatteringReaction(EqualityMixin):
r"""Thermal scattering reaction
This class is used to hold the integral and differential cross sections
for either elastic or inelastic thermal scattering.
Parameters
----------
xs : dict of str to Function1D
Integral cross section at each temperature
distribution : dict of str to AngleEnergy
Secondary angle-energy distribution at each temperature
Attributes
----------
xs : dict of str to Function1D
Integral cross section at each temperature
distribution : dict of str to AngleEnergy
Secondary angle-energy distribution at each temperature
"""
def __init__(self, xs, distribution):
self.xs = xs
self.distribution = distribution
def to_hdf5(self, group, name):
"""Write thermal scattering reaction to HDF5
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : {'elastic', 'inelastic'}
Name of reaction to write
"""
for T, xs in self.xs.items():
Tgroup = group.require_group(T)
rx_group = Tgroup.create_group(name)
xs.to_hdf5(rx_group, 'xs')
dgroup = rx_group.create_group('distribution')
self.distribution[T].to_hdf5(dgroup)
@classmethod
def from_hdf5(cls, group, name, temperatures):
"""Generate thermal scattering reaction data from HDF5
Parameters
----------
group : h5py.Group
HDF5 group to read from
name : {'elastic', 'inelastic'}
Name of the reaction to read
temperatures : Iterable of str
Temperatures to read
Returns
-------
openmc.data.ThermalScatteringReaction
Thermal scattering reaction data
"""
xs = {}
distribution = {}
for T in temperatures:
rx_group = group[T][name]
xs[T] = Function1D.from_hdf5(rx_group['xs'])
if isinstance(xs[T], CoherentElastic):
distribution[T] = CoherentElasticAE(xs[T])
else:
distribution[T] = AngleEnergy.from_hdf5(rx_group['distribution'])
return cls(xs, distribution)
class ThermalScattering(EqualityMixin):
"""A ThermalScattering object contains thermal scattering data as represented by
an S(alpha, beta) table.
@ -223,10 +394,11 @@ class ThermalScattering(EqualityMixin):
----------
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
elastic_xs : openmc.data.Tabulated1D or openmc.data.CoherentElastic
Elastic scattering cross section derived in the coherent or incoherent
approximation
inelastic_xs : openmc.data.Tabulated1D
energy_max : float
Maximum energy for thermal scattering data in [eV]
elastic : openmc.data.ThermalScatteringReaction or None
Elastic scattering derived in the coherent or incoherent approximation
inelastic : openmc.data.ThermalScatteringReaction
Inelastic scattering cross section derived in the incoherent
approximation
name : str
@ -243,28 +415,24 @@ class ThermalScattering(EqualityMixin):
"""
def __init__(self, name, atomic_weight_ratio, kTs):
def __init__(self, name, atomic_weight_ratio, energy_max, kTs):
self.name = name
self.atomic_weight_ratio = atomic_weight_ratio
self.energy_max = energy_max
self.kTs = kTs
self.elastic_xs = {}
self.elastic_mu_out = {}
self.inelastic_xs = {}
self.inelastic_e_out = {}
self.inelastic_mu_out = {}
self.inelastic_dist = {}
self.secondary_mode = None
self.elastic = None
self.inelastic = None
self.nuclides = []
def __repr__(self):
if hasattr(self, 'name'):
return "<Thermal Scattering Data: {0}>".format(self.name)
return "<Thermal Scattering Data: {}>".format(self.name)
else:
return "<Thermal Scattering Data>"
@property
def temperatures(self):
return ["{}K".format(int(round(kT / K_BOLTZMANN))) for kT in self.kTs]
return [_temperature_str(kT / K_BOLTZMANN) for kT in self.kTs]
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export table to an HDF5 file.
@ -282,43 +450,23 @@ class ThermalScattering(EqualityMixin):
"""
# Open file and write version
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_thermal')
f.attrs['version'] = np.array(HDF5_VERSION)
with h5py.File(str(path), mode, libver=libver) as f:
f.attrs['filetype'] = np.string_('data_thermal')
f.attrs['version'] = np.array(HDF5_VERSION)
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['nuclides'] = np.array(self.nuclides, dtype='S')
g.attrs['secondary_mode'] = np.string_(self.secondary_mode)
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
ktg.create_dataset(temperature, data=self.kTs[i])
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['energy_max'] = self.energy_max
g.attrs['nuclides'] = np.array(self.nuclides, dtype='S')
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
ktg.create_dataset(temperature, data=self.kTs[i])
for T in self.temperatures:
Tg = g.create_group(T)
# Write thermal elastic scattering
if self.elastic_xs:
elastic_group = Tg.create_group('elastic')
self.elastic_xs[T].to_hdf5(elastic_group, 'xs')
if self.elastic_mu_out:
elastic_group.create_dataset('mu_out',
data=self.elastic_mu_out[T])
# Write thermal inelastic scattering
if self.inelastic_xs:
inelastic_group = Tg.create_group('inelastic')
self.inelastic_xs[T].to_hdf5(inelastic_group, 'xs')
if self.secondary_mode in ('equal', 'skewed'):
inelastic_group.create_dataset('energy_out',
data=self.inelastic_e_out[T])
inelastic_group.create_dataset('mu_out',
data=self.inelastic_mu_out[T])
elif self.secondary_mode == 'continuous':
self.inelastic_dist[T].to_hdf5(inelastic_group)
f.close()
# Write elastic/inelastic reaction data
if self.elastic is not None:
self.elastic.to_hdf5(g, 'elastic')
self.inelastic.to_hdf5(g, 'inelastic')
def add_temperature_from_ace(self, ace_or_filename, name=None):
"""Add data to the ThermalScattering object from an ACE file at a
@ -356,20 +504,14 @@ class ThermalScattering(EqualityMixin):
self.kTs += data.kTs
# Add inelastic cross section and distributions
if strT in data.inelastic_xs:
self.inelastic_xs[strT] = data.inelastic_xs[strT]
if strT in data.inelastic_e_out:
self.inelastic_e_out[strT] = data.inelastic_e_out[strT]
if strT in data.inelastic_mu_out:
self.inelastic_mu_out[strT] = data.inelastic_mu_out[strT]
if strT in data.inelastic_dist:
self.inelastic_dist[strT] = data.inelastic_dist[strT]
if data.inelastic is not None:
self.inelastic.xs.update(data.inelastic.xs)
self.inelastic.distribution.update(data.inelastic.distribution)
# Add elastic cross sectoin and angular distribution
if strT in data.elastic_xs:
self.elastic_xs[strT] = data.elastic_xs[strT]
if strT in data.elastic_mu_out:
self.elastic_mu_out[strT] = data.elastic_mu_out[strT]
if data.elastic is not None:
self.elastic.xs.update(data.elastic.xs)
self.elastic.distribution.update(data.elastic.distribution)
@classmethod
def from_hdf5(cls, group_or_filename):
@ -411,46 +553,23 @@ class ThermalScattering(EqualityMixin):
name = group.name[1:]
atomic_weight_ratio = group.attrs['atomic_weight_ratio']
energy_max = group.attrs['energy_max']
kTg = group['kTs']
kTs = []
for temp in kTg:
kTs.append(kTg[temp][()])
temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs]
kTs = [dataset[()] for dataset in kTg.values()]
table = cls(name, atomic_weight_ratio, kTs)
table = cls(name, atomic_weight_ratio, energy_max, kTs)
table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']]
table.secondary_mode = group.attrs['secondary_mode'].decode()
# Read thermal elastic scattering
for T in temperatures:
Tgroup = group[T]
if 'elastic' in Tgroup:
elastic_group = Tgroup['elastic']
if 'elastic' in group[table.temperatures[0]]:
table.elastic = ThermalScatteringReaction.from_hdf5(
group, 'elastic', table.temperatures
)
# Cross section
elastic_xs_type = elastic_group['xs'].attrs['type'].decode()
if elastic_xs_type == 'Tabulated1D':
table.elastic_xs[T] = Tabulated1D.from_hdf5(
elastic_group['xs'])
elif elastic_xs_type == 'bragg':
table.elastic_xs[T] = CoherentElastic.from_hdf5(
elastic_group['xs'])
# Angular distribution
if 'mu_out' in elastic_group:
table.elastic_mu_out[T] = elastic_group['mu_out'][()]
# Read thermal inelastic scattering
if 'inelastic' in Tgroup:
inelastic_group = Tgroup['inelastic']
table.inelastic_xs[T] = Tabulated1D.from_hdf5(
inelastic_group['xs'])
if table.secondary_mode in ('equal', 'skewed'):
table.inelastic_e_out[T] = inelastic_group['energy_out'][()]
table.inelastic_mu_out[T] = inelastic_group['mu_out'][()]
elif table.secondary_mode == 'continuous':
table.inelastic_dist[T] = AngleEnergy.from_hdf5(
inelastic_group)
# Read thermal inelastic scattering
table.inelastic = ThermalScatteringReaction.from_hdf5(
group, 'inelastic', table.temperatures
)
return table
@ -485,41 +604,30 @@ class ThermalScattering(EqualityMixin):
# Assign temperature to the running list
kTs = [ace.temperature*EV_PER_MEV]
temperatures = [str(int(round(ace.temperature*EV_PER_MEV
/ K_BOLTZMANN))) + "K"]
table = cls(name, ace.atomic_weight_ratio, kTs)
# Incoherent inelastic scattering cross section
idx = ace.jxs[1]
n_energy = int(ace.xss[idx])
energy = ace.xss[idx+1 : idx+1+n_energy]*EV_PER_MEV
xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy]
table.inelastic_xs[temperatures[0]] = Tabulated1D(energy, xs)
if ace.nxs[7] == 0:
table.secondary_mode = 'equal'
elif ace.nxs[7] == 1:
table.secondary_mode = 'skewed'
elif ace.nxs[7] == 2:
table.secondary_mode = 'continuous'
inelastic_xs = Tabulated1D(energy, xs)
energy_max = energy[-1]
# Incoherent inelastic angle-energy distribution
continuous = (ace.nxs[7] == 2)
n_energy_out = ace.nxs[4]
if table.secondary_mode in ('equal', 'skewed'):
if not continuous:
n_mu = ace.nxs[3]
idx = ace.jxs[3]
table.inelastic_e_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2):
n_mu + 2]*EV_PER_MEV
table.inelastic_e_out[temperatures[0]].shape = \
(n_energy, n_energy_out)
energy_out = ace.xss[idx:idx + n_energy * n_energy_out *
(n_mu + 2): n_mu + 2]*EV_PER_MEV
energy_out.shape = (n_energy, n_energy_out)
table.inelastic_mu_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)]
table.inelastic_mu_out[temperatures[0]].shape = \
(n_energy, n_energy_out, n_mu+2)
table.inelastic_mu_out[temperatures[0]] = \
table.inelastic_mu_out[temperatures[0]][:, :, 1:]
mu_out = ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)]
mu_out.shape = (n_energy, n_energy_out, n_mu+2)
mu_out = mu_out[:, :, 1:]
skewed = (ace.nxs[7] == 1)
distribution = IncoherentInelasticAEDiscrete(energy_out, mu_out, skewed)
else:
n_mu = ace.nxs[3] - 1
idx = ace.jxs[3]
@ -558,10 +666,16 @@ class ThermalScattering(EqualityMixin):
# Create correlated angle-energy distribution
breakpoints = [n_energy]
interpolation = [2]
energy = table.inelastic_xs[temperatures[0]].x
table.inelastic_dist[temperatures[0]] = CorrelatedAngleEnergy(
energy = inelastic_xs.x
distribution = IncoherentInelasticAE(
breakpoints, interpolation, energy, energy_out, mu_out)
table = cls(name, ace.atomic_weight_ratio, energy_max, kTs)
T = table.temperatures[0]
table.inelastic = ThermalScatteringReaction(
{T: inelastic_xs}, {T: distribution}
)
# Incoherent/coherent elastic scattering cross section
idx = ace.jxs[4]
n_mu = ace.nxs[6] + 1
@ -572,22 +686,23 @@ class ThermalScattering(EqualityMixin):
if ace.nxs[5] == 4:
# Coherent elastic
table.elastic_xs[temperatures[0]] = CoherentElastic(
energy, P*EV_PER_MEV)
xs = CoherentElastic(energy, P*EV_PER_MEV)
distribution = CoherentElasticAE(xs)
# Coherent elastic shouldn't have angular distributions listed
assert n_mu == 0
else:
# Incoherent elastic
table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P)
xs = Tabulated1D(energy, P)
# Angular distribution
assert n_mu > 0
idx = ace.jxs[6]
table.elastic_mu_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_mu]
table.elastic_mu_out[temperatures[0]].shape = \
(n_energy, n_mu)
mu_out = ace.xss[idx:idx + n_energy * n_mu]
mu_out.shape = (n_energy, n_mu)
distribution = IncoherentElasticAEDiscrete(mu_out)
table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution})
# Get relevant nuclides -- NJOY only allows one to specify three
# nuclides that the S(a,b) table applies to. Thus, for all elements
@ -611,8 +726,9 @@ class ThermalScattering(EqualityMixin):
@classmethod
def from_njoy(cls, filename, filename_thermal, temperatures=None,
evaluation=None, evaluation_thermal=None, **kwargs):
"""Generate incident neutron data by running NJOY.
evaluation=None, evaluation_thermal=None,
use_endf_data=True, **kwargs):
"""Generate thermal scattering data by running NJOY.
Parameters
----------
@ -631,6 +747,9 @@ class ThermalScattering(EqualityMixin):
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to
use.
use_endf_data : bool
If the material has incoherent elastic scattering, the ENDF data
will be used rather than the ACE data.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal`
@ -655,4 +774,150 @@ class ThermalScattering(EqualityMixin):
for table in lib.tables[1:]:
data.add_temperature_from_ace(table)
# Load ENDF data to replace incoherent elastic
if use_endf_data:
data_endf = cls.from_endf(filename_thermal)
if data_endf.elastic is not None:
# Get appropriate temperatures
if temperatures is None:
temperatures = data_endf.temperatures
else:
temperatures = [_temperature_str(t) for t in temperatures]
# 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):
rx.xs[t] = rx_endf.xs[t]
rx.distribution[t] = rx_endf.distribution[t]
return data
@classmethod
def from_endf(cls, ev_or_filename):
"""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.
Returns
-------
openmc.data.ThermalScattering
Thermal scattering data
"""
if isinstance(ev_or_filename, endf.Evaluation):
ev = ev_or_filename
else:
ev = endf.Evaluation(ev_or_filename)
# 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
# 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])
elif lhtr == 2:
# incoherent elastic
params, W = endf.get_tab1_record(file_obj)
bound_xs = params[0]
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)
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])
params = endf.get_head_record(file_obj)
data = {'symmetric': params[4] == 0}
# Get information about principal atom
params, B = endf.get_list_record(file_obj)
data['log'] = bool(params[2])
data['free_atom_xs'] = B[0]
data['epsilon'] = B[1]
data['A0'] = awr = B[2]
data['e_max'] = energy_max = B[3]
data['M0'] = B[5]
# Get information about non-principal atoms
n_non_principal = params[5]
data['non_principal'] = []
NonPrincipal = namedtuple('NonPrincipal', ['func', 'xs', 'A', 'M'])
for i in range(1, n_non_principal + 1):
func = {0.0: 'SCT', 1.0: 'free gas', 2.0: 'diffusive'}[B[6*i]]
xs = B[6*i + 1]
A = B[6*i + 2]
M = B[6*i + 5]
data['non_principal'].append(NonPrincipal(func, xs, A, M))
# Get S(alpha,beta,T)
kTs = []
if data['free_atom_xs'] > 0.0:
params, _ = endf.get_tab2_record(file_obj)
n_beta = params[5]
sab = {'beta': np.empty(n_beta)}
for i in range(n_beta):
params, S = endf.get_tab1_record(file_obj)
T0, beta, lt = params[:3]
if i == 0:
sab['alpha'] = alpha = S.x
sab[T0] = np.empty((alpha.size, n_beta))
kTs.append(K_BOLTZMANN * T0)
sab['beta'][i] = beta
sab[T0][:, i] = S.y
for _ in range(lt):
params, S = endf.get_list_record(file_obj)
T = params[0]
if i == 0:
sab[T] = np.empty((alpha.size, n_beta))
kTs.append(K_BOLTZMANN * T)
sab[T][:, i] = S
data['sab'] = sab
# Get effective temperature for each atom
_, Teff = endf.get_tab1_record(file_obj)
data['effective_temperature'] = [Teff]
for atom in data['non_principal']:
if atom.func == 'SCT':
_, Teff = endf.get_tab1_record(file_obj)
data['effective_temperature'].append(Teff)
name = ev.target['zsymam'].strip()
instance = cls(name, awr, energy_max, kTs)
if elastic is not None:
instance.elastic = elastic
# Currently we don't have a proper cross section or distribution for
# incoherent inelastic, so we just create an empty object and attach
# all the data as a dictionary
instance.inelastic = ThermalScatteringReaction(None, None)
instance.inelastic.data = data
return instance

212
openmc/data/thermal_angle_energy.py Normal file
View file

@ -0,0 +1,212 @@
import numpy as np
from .angle_energy import AngleEnergy
from .correlated import CorrelatedAngleEnergy
class CoherentElasticAE(AngleEnergy):
r"""Differential cross section for coherent elastic scattering
The differential cross section for coherent elastic scattering from a
powdered crystalline material may be represented as:
.. math::
\frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{1}{E} \sum
\limits_{i=1}^{E_i < E} s_i(T) \delta(\mu - \mu_i) \delta (E - E')
/(2\pi)
where :math:`E_i` are the energies of the Bragg edges in [eV], :math:`s_i(T)`
is the structure factor in [eV-b] at the moderator temperature :math:`T`
in [K], and :math:`\mu_i = 1 - 2E_i/E`.
Parameters
----------
coherent_xs : openmc.data.CoherentElastic
Coherent elastic scattering cross section
Attributes
----------
coherent_xs : openmc.data.CoherentElastic
Coherent elastic scattering cross section
"""
def __init__(self, coherent_xs):
self.coherent_xs = coherent_xs
def to_hdf5(self, group):
"""Write coherent elastic distribution to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['type'] = np.string_('coherent_elastic')
group['coherent_xs'] = group.parent['xs']
class IncoherentElasticAE(AngleEnergy):
r"""Differential cross section for incoherent elastic scattering
The differential cross section for incoherent elastic scattering may be
represented as:
.. math::
\frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{\sigma_b}
{4\pi} e^{-2EW'(T)(1-\mu)} \delta(E - E')
where :math:`\sigma_b` is the characteristic cross section in [b] and
:math:`W'(T)` is the Debye-Waller integral divided by the atomic mass in
[eV\ :math:`^{-1}`].
Parameters
----------
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
Attributes
----------
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
"""
def __init__(self, debye_waller):
self.debye_waller = debye_waller
def to_hdf5(self, group):
"""Write incoherent elastic distribution to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['type'] = np.string_('incoherent_elastic')
group.create_dataset('debye_waller', data=self.debye_waller)
@classmethod
def from_hdf5(cls, group):
"""Generate incoherent elastic distribution from HDF5 data
Parameters
----------
group : h5py.Group
HDF5 group to read from
Returns
-------
openmc.data.IncoherentElasticAE
Incoherent elastic distribution
"""
return cls(group['debye_waller'])
class IncoherentElasticAEDiscrete(AngleEnergy):
"""Discrete angle representation of incoherent elastic scattering
Parameters
----------
mu_out : numpy.ndarray
Equi-probable discrete angles at each incoming energy
"""
def __init__(self, mu_out):
self.mu_out = mu_out
def to_hdf5(self, group):
"""Write discrete incoherent elastic distribution to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['type'] = np.string_('incoherent_elastic_discrete')
group.create_dataset('mu_out', data=self.mu_out)
@classmethod
def from_hdf5(cls, group):
"""Generate discrete incoherent elastic distribution from HDF5 data
Parameters
----------
group : h5py.Group
HDF5 group to read from
Returns
-------
openmc.data.IncoherentElasticAEDiscrete
Discrete incoherent elastic distribution
"""
return cls(group['mu_out'][()])
class IncoherentInelasticAEDiscrete(AngleEnergy):
"""Discrete angle representation of incoherent inelastic scattering
Parameters
----------
energy_out : numpy.ndarray
Outgoing energies for each incoming energy
mu_out : numpy.ndarray
Discrete angles for each incoming/outgoing energy
skewed : bool
Whether discrete angles are equi-probable or have a skewed distribution
Attributes
----------
energy_out : numpy.ndarray
Outgoing energies for each incoming energy
mu_out : numpy.ndarray
Discrete angles for each incoming/outgoing energy
skewed : bool
Whether discrete angles are equi-probable or have a skewed distribution
"""
def __init__(self, energy_out, mu_out, skewed=False):
self.energy_out = energy_out
self.mu_out = mu_out
self.skewed = skewed
def to_hdf5(self, group):
"""Write discrete incoherent inelastic distribution to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['type'] = np.string_('incoherent_inelastic_discrete')
group.create_dataset('energy_out', data=self.energy_out)
group.create_dataset('mu_out', data=self.mu_out)
group.create_dataset('skewed', data=self.skewed)
@classmethod
def from_hdf5(cls, group):
"""Generate discrete incoherent inelastic distribution from HDF5 data
Parameters
----------
group : h5py.Group
HDF5 group to read from
Returns
-------
openmc.data.IncoherentInelasticAEDiscrete
Discrete incoherent inelastic distribution
"""
energy_out = group['energy_out'][()]
mu_out = group['mu_out'][()]
skewed = bool(group['skewed'])
return cls(energy_out, mu_out, skewed)
class IncoherentInelasticAE(CorrelatedAngleEnergy):
_name = 'incoherent_inelastic'

7
openmc/deplete/abc.py
View file

@ -49,6 +49,9 @@ class TransportOperator(metaclass=ABCMeta):
Path to the depletion chain XML file. Defaults to the file
listed under ``depletion_chain`` in
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``.
Attributes
----------
@ -58,7 +61,7 @@ class TransportOperator(metaclass=ABCMeta):
nuclides with reaction rates. Defaults to 1.0e3.
"""
def __init__(self, chain_file=None):
def __init__(self, chain_file=None, fission_q=None):
self.dilute_initial = 1.0e3
self.output_dir = '.'
@ -81,7 +84,7 @@ class TransportOperator(metaclass=ABCMeta):
warn("Use of OPENMC_DEPLETE_CHAIN is deprecated in favor "
"of adding depletion_chain to OPENMC_CROSS_SECTIONS",
FutureWarning)
self.chain = Chain.from_xml(chain_file)
self.chain = Chain.from_xml(chain_file, fission_q)
@abstractmethod
def __call__(self, vec, print_out=True):

24
openmc/deplete/chain.py
View file

@ -4,16 +4,18 @@ This module contains information about a depletion chain. A depletion chain is
loaded from an .xml file and all the nuclides are linked together.
"""
from collections import OrderedDict, defaultdict
from io import StringIO
from itertools import chain
import math
import re
import os
from collections import OrderedDict, defaultdict
from collections.abc import Mapping
from openmc.checkvalue import check_type
# Try to use lxml if it is available. It preserves the order of attributes and
# provides a pretty-printer by default. If not available, use OpenMC function to
# pretty print.
# provides a pretty-printer by default. If not available,
# use OpenMC function to pretty print.
try:
import lxml.etree as ET
_have_lxml = True
@ -312,22 +314,32 @@ class Chain(object):
return chain
@classmethod
def from_xml(cls, filename):
def from_xml(cls, filename, fission_q=None):
"""Reads a depletion chain XML file.
Parameters
----------
filename : str
The path to the depletion chain XML file.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
If not given, values will be pulled from ``filename``
"""
chain = cls()
if fission_q is not None:
check_type("fission_q", fission_q, Mapping)
else:
fission_q = {}
# Load XML tree
root = ET.parse(str(filename))
for i, nuclide_elem in enumerate(root.findall('nuclide')):
nuc = Nuclide.from_xml(nuclide_elem)
this_q = fission_q.get(nuclide_elem.get("name"))
nuc = Nuclide.from_xml(nuclide_elem, this_q)
chain.nuclide_dict[nuc.name] = i
# Check for reaction paths

7
openmc/deplete/nuclide.py
View file

@ -112,13 +112,16 @@ class Nuclide(object):
return len(self.reactions)
@classmethod
def from_xml(cls, element):
def from_xml(cls, element, fission_q=None):
"""Read nuclide from an XML element.
Parameters
----------
element : xml.etree.ElementTree.Element
XML element to write nuclide data to
fission_q : None or float
User-supplied fission Q value [eV].
Will be read from the element if not given
Returns
-------
@ -153,6 +156,8 @@ class Nuclide(object):
target = reaction_elem.get('target')
else:
target = None
if fission_q is not None:
Q = fission_q
# Append reaction
nuc.reactions.append(ReactionTuple(

10
openmc/deplete/operator.py
View file

@ -73,6 +73,9 @@ class Operator(TransportOperator):
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV]. If not given,
values will be pulled from the ``chain_file``.
Attributes
----------
@ -110,14 +113,14 @@ class Operator(TransportOperator):
"""
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False):
super().__init__(chain_file)
diff_burnable_mats=False, fission_q=None):
super().__init__(chain_file, fission_q)
self.round_number = False
self.settings = settings
self.geometry = geometry
self.diff_burnable_mats = diff_burnable_mats
if prev_results != None:
if prev_results is not None:
# Reload volumes into geometry
prev_results[-1].transfer_volumes(geometry)
@ -149,7 +152,6 @@ class Operator(TransportOperator):
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
def __call__(self, vec, power, print_out=True):
"""Runs a simulation.

10
openmc/model/model.py
View file

@ -124,7 +124,7 @@ class Model(object):
self._plots.append(plot)
def deplete(self, timesteps, chain_file=None, method='cecm',
**kwargs):
fission_q=None, **kwargs):
"""Deplete model using specified timesteps/power
Parameters
@ -138,6 +138,9 @@ class Model(object):
:envvar:`OPENMC_CROSS_SECTIONS` environment variable if it exists.
method : str
Integration method used for depletion (e.g., 'cecm', 'predictor')
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
If not given, values will be pulled from the ``chain_file``.
**kwargs
Keyword arguments passed to integration function (e.g.,
:func:`openmc.deplete.integrator.cecm`)
@ -149,7 +152,10 @@ class Model(object):
import openmc.deplete as dep
# Create OpenMC transport operator
op = dep.Operator(self.geometry, self.settings, chain_file)
op = dep.Operator(
self.geometry, self.settings, chain_file,
fission_q=fission_q,
)
# Perform depletion
check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4',

24
src/endf.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/endf.h"
#include <algorithm> // for copy
#include <array>
#include <cmath> // for log, exp
#include <iterator> // for back_inserter
#include <stdexcept> // for runtime_error
@ -87,6 +88,10 @@ read_function(hid_t group, const char* name)
func = std::make_unique<Tabulated1D>(dset);
} else if (func_type == "Polynomial") {
func = std::make_unique<Polynomial>(dset);
} else if (func_type == "CoherentElastic") {
func = std::make_unique<CoherentElasticXS>(dset);
} else if (func_type == "IncoherentElastic") {
func = std::make_unique<IncoherentElasticXS>(dset);
} else {
throw std::runtime_error{"Unknown function type " + func_type +
" for dataset " + object_name(dset)};
@ -234,4 +239,23 @@ double CoherentElasticXS::operator()(double E) const
}
}
//==============================================================================
// IncoherentElasticXS implementation
//==============================================================================
IncoherentElasticXS::IncoherentElasticXS(hid_t dset)
{
std::array<double, 2> tmp;
read_dataset(dset, nullptr, tmp);
bound_xs_ = tmp[0];
debye_waller_ = tmp[1];
}
double IncoherentElasticXS::operator()(double E) const
{
// Determine cross section using ENDF-102, Eq. (7.5)
double W = debye_waller_;
return bound_xs_ / 2.0 * ((1 - std::exp(-4.0*E*W))/(2.0*E*W));
}
} // namespace openmc

2
src/material.cpp
View file

@ -774,7 +774,7 @@ void Material::calculate_neutron_xs(Particle& p) const
// If particle energy is greater than the highest energy for the
// S(a,b) table, then don't use the S(a,b) table
if (p.E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
if (p.E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE;
// Increment position in thermal_tables_
++j;

2
src/reaction.cpp
View file

@ -42,8 +42,6 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
// Get threshold index
TemperatureXS xs;
read_attribute(dset, "threshold_idx", xs.threshold);
// TODO: change HDF5 format so that threshold_idx is 0-based
--xs.threshold;
// Read cross section values
read_dataset(dset, xs.value);

10
src/reaction_product.cpp
View file

@ -3,6 +3,7 @@
#include <memory> // for unique_ptr
#include <string> // for string
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/secondary_correlated.h"
@ -42,14 +43,7 @@ ReactionProduct::ReactionProduct(hid_t group)
read_attribute(group, "decay_rate", decay_rate_);
// Read secondary particle yield
hid_t yield = open_dataset(group, "yield");
read_attribute(yield, "type", temp);
if (temp == "Tabulated1D") {
yield_ = std::unique_ptr<Function1D>{new Tabulated1D{yield}};
} else if (temp == "Polynomial") {
yield_ = std::unique_ptr<Function1D>{new Polynomial{yield}};
}
close_dataset(yield);
yield_ = read_function(group, "yield");
int n;
read_attribute(group, "n_distribution", n);

327
src/secondary_thermal.cpp Normal file
View file

@ -0,0 +1,327 @@
#include "openmc/secondary_thermal.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "xtensor/xview.hpp"
#include <cmath> // for log, exp
namespace openmc {
// Helper function to get index on incident energy grid
void
get_energy_index(const std::vector<double>& energies, double E, int& i, double& f)
{
// Get index and interpolation factor for elastic grid
i = 0;
f = 0.0;
if (E >= energies.front()) {
i = lower_bound_index(energies.begin(), energies.end(), E);
f = (E - energies[i]) / (energies[i+1] - energies[i]);
}
}
//==============================================================================
// CoherentElasticAE implementation
//==============================================================================
CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs)
: xs_{xs}
{ }
void
CoherentElasticAE::sample(double E_in, double& E_out, double& mu) const
{
// Get index and interpolation factor for elastic grid
int i;
double f;
const auto& energies {xs_.bragg_edges()};
get_energy_index(energies, E_in, i, f);
// Sample a Bragg edge between 1 and i
const auto& factors = xs_.factors();
double prob = prn() * factors[i+1];
int k = 0;
if (prob >= factors.front()) {
k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob);
}
// Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
mu = 1.0 - 2.0*energies[k] / E_in;
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
E_out = E_in;
}
//==============================================================================
// IncoherentElasticAE implementation
//==============================================================================
IncoherentElasticAE::IncoherentElasticAE(hid_t group)
{
read_attribute(group, "debye_waller", debye_waller_);
}
void
IncoherentElasticAE::sample(double E_in, double& E_out, double& mu) const
{
// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
double c = 2 * E_in * debye_waller_;
mu = std::log(1.0 + prn()*(std::exp(2.0*c) - 1))/c - 1.0;
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
E_out = E_in;
}
//==============================================================================
// IncoherentElasticAEDiscrete implementation
//==============================================================================
IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
const std::vector<double>& energy)
: energy_{energy}
{
read_dataset(group, "mu_out", mu_out_);
}
void
IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const
{
// Get index and interpolation factor for elastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
// Interpolate between two discrete cosines corresponding to neighboring
// incoming energies.
// Sample outgoing cosine bin
int k = prn() * mu_out_.shape()[1];
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ik = mu_out_(i, k);
double mu_i1k = mu_out_(i+1, k);
// Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ik + f*mu_i1k;
// Energy doesn't change in elastic scattering
E_out = E_in;
}
//==============================================================================
// IncoherentInelasticAEDiscrete implementation
//==============================================================================
IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
const std::vector<double>& energy)
: energy_{energy}
{
read_dataset(group, "energy_out", energy_out_);
read_dataset(group, "mu_out", mu_out_);
read_dataset(group, "skewed", skewed_);
}
void
IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const
{
// Get index and interpolation factor for inelastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
// Now that we have an incoming energy bin, we need to determine the outgoing
// energy bin. This will depend on whether the outgoing energy distribution is
// skewed. If it is skewed, then the first two and last two bins have lower
// probabilities than the other bins (0.1 for the first and last bins and 0.4
// for the second and second to last bins, relative to a normal bin
// probability of 1). Otherwise, each bin is equally probable.
int j;
int n = energy_out_.shape()[1];
if (!skewed_) {
// All bins equally likely
j = prn() * n;
} else {
// Distribution skewed away from edge points
double r = prn() * (n - 3);
if (r > 1.0) {
// equally likely N-4 middle bins
j = r + 1;
} else if (r > 0.6) {
// second to last bin has relative probability of 0.4
j = n - 2;
} else if (r > 0.5) {
// last bin has relative probability of 0.1
j = n - 1;
} else if (r > 0.1) {
// second bin has relative probability of 0.4
j = 1;
} else {
// first bin has relative probability of 0.1
j = 0;
}
}
// Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
double E_ij = energy_out_(i, j);
double E_i1j = energy_out_(i+1, j);
// Outgoing energy
E_out = (1 - f)*E_ij + f*E_i1j;
// Sample outgoing cosine bin
int m = mu_out_.shape()[2];
int k = prn() * m;
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ijk = mu_out_(i, j, k);
double mu_i1jk = mu_out_(i+1, j, k);
// Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk;
}
//==============================================================================
// IncoherentInelasticAE implementation
//==============================================================================
IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
{
// Read correlated angle-energy distribution
CorrelatedAngleEnergy dist {group};
// Copy incident energies
energy_ = dist.energy();
// Convert to S(a,b) native format
for (const auto& edist : dist.distribution()) {
// Create temporary distribution
DistEnergySab d;
// Copy outgoing energy distribution
d.n_e_out = edist.e_out.size();
d.e_out = edist.e_out;
d.e_out_pdf = edist.p;
d.e_out_cdf = edist.c;
for (int j = 0; j < d.n_e_out; ++j) {
auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());
if (adist) {
// On first pass, allocate space for angles
if (j == 0) {
auto n_mu = adist->x().size();
d.mu = xt::empty<double>({d.n_e_out, n_mu});
}
// Copy outgoing angles
auto mu_j = xt::view(d.mu, j);
std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());
}
}
distribution_.emplace_back(std::move(d));
}
}
void
IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const
{
// Get index and interpolation factor for inelastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
// Sample between ith and [i+1]th bin
int l = f > prn() ? i + 1 : i;
// Determine endpoints on grid i
auto n = distribution_[i].e_out.size();
double E_i_1 = distribution_[i].e_out[0];
double E_i_J = distribution_[i].e_out[n - 1];
// Determine endpoints on grid i + 1
n = distribution_[i + 1].e_out.size();
double E_i1_1 = distribution_[i + 1].e_out[0];
double E_i1_J = distribution_[i + 1].e_out[n - 1];
double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
double E_J = E_i_J + f * (E_i1_J - E_i_J);
// Determine outgoing energy bin
// (First reset n_energy_out to the right value)
n = distribution_[l].n_e_out;
double r1 = prn();
double c_j = distribution_[l].e_out_cdf[0];
double c_j1;
std::size_t j;
for (j = 0; j < n - 1; ++j) {
c_j1 = distribution_[l].e_out_cdf[j + 1];
if (r1 < c_j1) break;
c_j = c_j1;
}
// check to make sure j is <= n_energy_out - 2
j = std::min(j, n - 2);
// Get the data to interpolate between
double E_l_j = distribution_[l].e_out[j];
double p_l_j = distribution_[l].e_out_pdf[j];
// Next part assumes linear-linear interpolation in standard
double E_l_j1 = distribution_[l].e_out[j + 1];
double p_l_j1 = distribution_[l].e_out_pdf[j + 1];
// Find secondary energy (variable E)
double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
if (frac == 0.0) {
E_out = E_l_j + (r1 - c_j) / p_l_j;
} else {
E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j +
2.0*frac*(r1 - c_j))) - p_l_j) / frac;
}
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
} else {
E_out = E_1 + (E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
}
// Sample outgoing cosine bin
int n_mu = distribution_[l].mu.shape()[1];
std::size_t k = prn() * n_mu;
// Rather than use the sampled discrete mu directly, it is smeared over
// a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// discrete mu value itself.
const auto& mu_l = distribution_[l].mu;
f = (r1 - c_j)/(c_j1 - c_j);
// Determine (k-1)th mu value
double mu_left;
if (k == 0) {
mu_left = -1.0;
} else {
mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
}
// Determine kth mu value
mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k));
// Determine (k+1)th mu value
double mu_right;
if (k == n_mu - 1) {
mu_right = 1.0;
} else {
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
}
// Smear angle
mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
}
} // namespace openmc

417
src/thermal.cpp
View file

@ -12,10 +12,12 @@
#include "xtensor/xview.hpp"
#include "openmc/constants.h"
#include "openmc/endf.h"
#include "openmc/error.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/secondary_correlated.h"
#include "openmc/secondary_thermal.h"
#include "openmc/settings.h"
namespace openmc {
@ -42,17 +44,8 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
name_ = name_.substr(1);
read_attribute(group, "atomic_weight_ratio", awr_);
read_attribute(group, "energy_max", energy_max_);
read_attribute(group, "nuclides", nuclides_);
std::string sec_mode;
read_attribute(group, "secondary_mode", sec_mode);
int secondary_mode;
if (sec_mode == "equal") {
secondary_mode = SAB_SECONDARY_EQUAL;
} else if (sec_mode == "skewed") {
secondary_mode = SAB_SECONDARY_SKEWED;
} else if (sec_mode == "continuous") {
secondary_mode = SAB_SECONDARY_CONT;
}
// Read temperatures
hid_t kT_group = open_group(group, "kTs");
@ -148,7 +141,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
// Open group for temperature i
hid_t T_group = open_group(group, temp_str.data());
data_.emplace_back(T_group, secondary_mode);
data_.emplace_back(T_group);
close_group(T_group);
}
@ -185,58 +178,8 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
// Set temperature index
*i_temp = i;
// Get pointer to S(a,b) table
auto& sab = data_[i];
// Get index and interpolation factor for inelastic grid
int i_grid = 0;
double f = 0.0;
if (E >= sab.inelastic_e_in_.front()) {
auto& E_in = sab.inelastic_e_in_;
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]);
}
// Calculate S(a,b) inelastic scattering cross section
auto& xs = sab.inelastic_sigma_;
*inelastic = xs[i_grid] + f * (xs[i_grid + 1] - xs[i_grid]);
// Check for elastic data
if (!sab.elastic_e_in_.empty()) {
// Determine whether elastic scattering is given in the coherent or
// incoherent approximation. For coherent, the cross section is
// represented as P/E whereas for incoherent, it is simply P
auto& E_in = sab.elastic_e_in_;
if (sab.elastic_mode_ == SAB_ELASTIC_COHERENT) {
if (E < E_in.front()) {
// If energy is below that of the lowest Bragg peak, the elastic
// cross section will be zero
*elastic = 0.0;
} else {
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
*elastic = sab.elastic_P_[i_grid] / E;
}
} else {
// Determine index on elastic energy grid
if (E < E_in.front()) {
i_grid = 0;
} else {
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
}
// Get interpolation factor for elastic grid
f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]);
// Calculate S(a,b) elastic scattering cross section
auto& xs = sab.elastic_P_;
*elastic = xs[i_grid] + f*(xs[i_grid + 1] - xs[i_grid]);
}
} else {
// No elastic data
*elastic = 0.0;
}
// Calculate cross sections for ith temperature
data_[i].calculate_xs(E, elastic, inelastic);
}
bool
@ -250,46 +193,39 @@ ThermalScattering::has_nuclide(const char* name) const
// ThermalData implementation
//==============================================================================
ThermalData::ThermalData(hid_t group, int secondary_mode)
: inelastic_mode_{secondary_mode}
ThermalData::ThermalData(hid_t group)
{
// Coherent elastic data
// Coherent/incoherent elastic data
if (object_exists(group, "elastic")) {
// Read cross section data
hid_t elastic_group = open_group(group, "elastic");
// Read elastic cross section
xt::xarray<double> temp;
hid_t dset = open_dataset(elastic_group, "xs");
read_dataset(dset, temp);
elastic_.xs = read_function(elastic_group, "xs");
// Get view on energies and cross section/probability values
auto E_in = xt::view(temp, 0);
auto P = xt::view(temp, 1);
// Read angle-energy distribution
hid_t dgroup = open_group(elastic_group, "distribution");
std::string temp;
read_attribute(dgroup, "type", temp);
if (temp == "coherent_elastic") {
auto xs = dynamic_cast<CoherentElasticXS*>(elastic_.xs.get());
elastic_.distribution = std::make_unique<CoherentElasticAE>(*xs);
// Set cross section data and type
std::copy(E_in.begin(), E_in.end(), std::back_inserter(elastic_e_in_));
std::copy(P.begin(), P.end(), std::back_inserter(elastic_P_));
n_elastic_e_in_ = elastic_e_in_.size();
// Set threshold energy
threshold_elastic_ = xs->bragg_edges().back();
// Determine whether elastic scattering is incoherent or coherent
std::string type;
read_attribute(dset, "type", type);
if (type == "tab1") {
elastic_mode_ = SAB_ELASTIC_INCOHERENT;
} else if (type == "bragg") {
elastic_mode_ = SAB_ELASTIC_COHERENT;
}
close_dataset(dset);
} else {
auto xs = dynamic_cast<Tabulated1D*>(elastic_.xs.get());
if (temp == "incoherent_elastic") {
elastic_.distribution = std::make_unique<IncoherentElasticAE>(dgroup);
} else if (temp == "incoherent_elastic_discrete") {
elastic_.distribution = std::make_unique<IncoherentElasticAEDiscrete>(
dgroup, xs->x()
);
}
// Set elastic threshold
threshold_elastic_ = elastic_e_in_.back();
// Read angle distribution
if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) {
xt::xarray<double> mu_out;
read_dataset(elastic_group, "mu_out", mu_out);
elastic_mu_ = mu_out;
// Set threshold energy
threshold_elastic_ = xs->x().back();
}
close_group(elastic_group);
@ -300,291 +236,58 @@ ThermalData::ThermalData(hid_t group, int secondary_mode)
// Read type of inelastic data
hid_t inelastic_group = open_group(group, "inelastic");
// Read cross section data
xt::xarray<double> temp;
read_dataset(inelastic_group, "xs", temp);
// Get view of inelastic cross section and energy grid
auto E_in = xt::view(temp, 0);
auto xs = xt::view(temp, 1);
// Set cross section data
std::copy(E_in.begin(), E_in.end(), std::back_inserter(inelastic_e_in_));
std::copy(xs.begin(), xs.end(), std::back_inserter(inelastic_sigma_));
n_inelastic_e_in_ = inelastic_e_in_.size();
// Read inelastic cross section
inelastic_.xs = read_function(inelastic_group, "xs");
// Set inelastic threshold
threshold_inelastic_ = inelastic_e_in_.back();
auto xs = dynamic_cast<Tabulated1D*>(inelastic_.xs.get());
threshold_inelastic_ = xs->x().back();
if (secondary_mode != SAB_SECONDARY_CONT) {
// Read energy distribution
xt::xarray<double> E_out;
read_dataset(inelastic_group, "energy_out", E_out);
inelastic_e_out_ = E_out;
n_inelastic_e_out_ = inelastic_e_out_.shape()[1];
// Read angle distribution
xt::xarray<double> mu_out;
read_dataset(inelastic_group, "mu_out", mu_out);
inelastic_mu_ = mu_out;
n_inelastic_mu_ = inelastic_mu_.shape()[2];
} else {
// Read correlated angle-energy distribution
CorrelatedAngleEnergy dist {inelastic_group};
// Convert to S(a,b) native format
for (const auto& edist : dist.distribution()) {
// Create temporary distribution
DistEnergySab d;
// Copy outgoing energy distribution
d.n_e_out = edist.e_out.size();
d.e_out = edist.e_out;
d.e_out_pdf = edist.p;
d.e_out_cdf = edist.c;
for (int j = 0; j < d.n_e_out; ++j) {
auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());
if (adist) {
// On first pass, allocate space for angles
if (j == 0) {
auto n_mu = adist->x().size();
n_inelastic_mu_ = n_mu;
d.mu = xt::empty<double>({d.n_e_out, n_mu});
}
// Copy outgoing angles
auto mu_j = xt::view(d.mu, j);
std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());
}
}
inelastic_data_.push_back(std::move(d));
}
// Read angle-energy distribution
hid_t dgroup = open_group(inelastic_group, "distribution");
std::string temp;
read_attribute(dgroup, "type", temp);
if (temp == "incoherent_inelastic") {
inelastic_.distribution = std::make_unique<IncoherentInelasticAE>(dgroup);
} else if (temp == "incoherent_inelastic_discrete") {
inelastic_.distribution = std::make_unique<IncoherentInelasticAEDiscrete>(
dgroup, xs->x()
);
}
close_group(inelastic_group);
}
}
void
ThermalData::calculate_xs(double E, double* elastic, double* inelastic) const
{
// Calculate thermal elastic scattering cross section
if (elastic_.xs) {
*elastic = (*elastic_.xs)(E);
} else {
*elastic = 0.0;
}
// Calculate thermal inelastic scattering cross section
*inelastic = (*inelastic_.xs)(E);
}
void
ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
double* E_out, double* mu)
{
// Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs.thermal_elastic / micro_xs.thermal) {
// elastic scattering
// Get index and interpolation factor for elastic grid
int i = 0;
double f = 0.0;
if (E >= elastic_e_in_.front()) {
auto& E_in = elastic_e_in_;
i = lower_bound_index(E_in.begin(), E_in.end(), E);
f = (E - E_in[i]) / (E_in[i+1] - E_in[i]);
}
// Select treatment based on elastic mode
if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) {
// With this treatment, we interpolate between two discrete cosines
// corresponding to neighboring incoming energies. This is used for
// data derived in the incoherent approximation
// Sample outgoing cosine bin
int k = prn() * n_elastic_mu_;
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ijk = elastic_mu_(i, k);
double mu_i1jk = elastic_mu_(i+1, k);
// Cosine of angle between incoming and outgoing neutron
*mu = (1 - f)*mu_ijk + f*mu_i1jk;
} else if (elastic_mode_ == SAB_ELASTIC_COHERENT) {
// This treatment is used for data derived in the coherent
// approximation, i.e. for crystalline structures that have Bragg
// edges.
// Sample a Bragg edge between 1 and i
double prob = prn() * elastic_P_[i+1];
int k = 0;
if (prob >= elastic_P_.front()) {
k = lower_bound_index(elastic_P_.begin(), elastic_P_.begin() + (i+1), prob);
}
// Characteristic scattering cosine for this Bragg edge
*mu = 1.0 - 2.0*elastic_e_in_[k] / E;
}
// Outgoing energy is same as incoming energy
*E_out = E;
elastic_.distribution->sample(E, *E_out, *mu);
} else {
// Perform inelastic calculations
// Get index and interpolation factor for inelastic grid
int i = 0;
double f = 0.0;
if (E >= inelastic_e_in_.front()) {
auto& E_in = inelastic_e_in_;
i = lower_bound_index(E_in.begin(), E_in.end(), E);
f = (E - E_in[i]) / (E_in[i+1] - E_in[i]);
}
// Now that we have an incoming energy bin, we need to determine the
// outgoing energy bin. This will depend on the "secondary energy
// mode". If the mode is 0, then the outgoing energy bin is chosen from a
// set of equally-likely bins. If the mode is 1, then the first
// two and last two bins are skewed to have lower probabilities than the
// other bins (0.1 for the first and last bins and 0.4 for the second and
// second to last bins, relative to a normal bin probability of 1).
// Finally, if the mode is 2, then a continuous distribution (with
// accompanying PDF and CDF is utilized)
if (inelastic_mode_ == SAB_SECONDARY_EQUAL ||
inelastic_mode_ == SAB_SECONDARY_SKEWED) {
int j;
if (inelastic_mode_ == SAB_SECONDARY_EQUAL) {
// All bins equally likely
j = prn() * n_inelastic_e_out_;
} else if (inelastic_mode_ == SAB_SECONDARY_SKEWED) {
// Distribution skewed away from edge points
double r = prn() * (n_inelastic_e_out_ - 3);
if (r > 1.0) {
// equally likely N-4 middle bins
j = r + 1;
} else if (r > 0.6) {
// second to last bin has relative probability of 0.4
j = n_inelastic_e_out_ - 2;
} else if (r > 0.5) {
// last bin has relative probability of 0.1
j = n_inelastic_e_out_ - 1;
} else if (r > 0.1) {
// second bin has relative probability of 0.4
j = 1;
} else {
// first bin has relative probability of 0.1
j = 0;
}
}
// Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
double E_ij = inelastic_e_out_(i, j);
double E_i1j = inelastic_e_out_(i+1, j);
// Outgoing energy
*E_out = (1 - f)*E_ij + f*E_i1j;
// Sample outgoing cosine bin
int k = prn() * n_inelastic_mu_;
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ijk = inelastic_mu_(i, j, k);
double mu_i1jk = inelastic_mu_(i+1, j, k);
// Cosine of angle between incoming and outgoing neutron
*mu = (1 - f)*mu_ijk + f*mu_i1jk;
} else if (inelastic_mode_ == SAB_SECONDARY_CONT) {
// Continuous secondary energy - this is to be similar to
// Law 61 interpolation on outgoing energy
// Sample between ith and [i+1]th bin
int l = f > prn() ? i + 1 : i;
// Determine endpoints on grid i
auto n = inelastic_data_[i].e_out.size();
double E_i_1 = inelastic_data_[i].e_out(0);
double E_i_J = inelastic_data_[i].e_out(n - 1);
// Determine endpoints on grid i + 1
n = inelastic_data_[i + 1].e_out.size();
double E_i1_1 = inelastic_data_[i + 1].e_out(0);
double E_i1_J = inelastic_data_[i + 1].e_out(n - 1);
double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
double E_J = E_i_J + f * (E_i1_J - E_i_J);
// Determine outgoing energy bin
// (First reset n_energy_out to the right value)
n = inelastic_data_[l].n_e_out;
double r1 = prn();
double c_j = inelastic_data_[l].e_out_cdf[0];
double c_j1;
std::size_t j;
for (j = 0; j < n - 1; ++j) {
c_j1 = inelastic_data_[l].e_out_cdf[j + 1];
if (r1 < c_j1) break;
c_j = c_j1;
}
// check to make sure j is <= n_energy_out - 2
j = std::min(j, n - 2);
// Get the data to interpolate between
double E_l_j = inelastic_data_[l].e_out[j];
double p_l_j = inelastic_data_[l].e_out_pdf[j];
// Next part assumes linear-linear interpolation in standard
double E_l_j1 = inelastic_data_[l].e_out[j + 1];
double p_l_j1 = inelastic_data_[l].e_out_pdf[j + 1];
// Find secondary energy (variable E)
double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
if (frac == 0.0) {
*E_out = E_l_j + (r1 - c_j) / p_l_j;
} else {
*E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j +
2.0*frac*(r1 - c_j))) - p_l_j) / frac;
}
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
*E_out = E_1 + (*E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
} else {
*E_out = E_1 + (*E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
}
// Sample outgoing cosine bin
std::size_t k = prn() * n_inelastic_mu_;
// Rather than use the sampled discrete mu directly, it is smeared over
// a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// discrete mu value itself.
const auto& mu_l = inelastic_data_[l].mu;
f = (r1 - c_j)/(c_j1 - c_j);
// Determine (k-1)th mu value
double mu_left;
if (k == 0) {
mu_left = -1.0;
} else {
mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
}
// Determine kth mu value
*mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k));
// Determine (k+1)th mu value
double mu_right;
if (k == n_inelastic_mu_ - 1) {
mu_right = 1.0;
} else {
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1));
}
// Smear angle
*mu += std::min(*mu - mu_left, mu_right - *mu)*(prn() - 0.5);
} // (inelastic secondary energy treatment)
} // (elastic or inelastic)
inelastic_.distribution->sample(E, *E_out, *mu);
}
// Because of floating-point roundoff, it may be possible for mu to be
// outside of the range [-1,1). In these cases, we just set mu to exactly
// -1 or 1
if (std::abs(*mu) > 1.0) *mu = std::copysign(1.0, *mu);
}
void free_memory_thermal()

179
tests/unit_tests/test_data_thermal.py
View file

@ -1,7 +1,9 @@
#!/usr/bin/env python
from collections.abc import Callable
from math import exp
import os
import random
import numpy as np
import pytest
@ -29,23 +31,53 @@ def graphite():
@pytest.fixture(scope='module')
def h2o_njoy():
"""H in H2O generated using NJOY."""
path_h1 = os.path.join(_ENDF_DATA, 'neutrons', 'n-001_H_001.endf')
path_h2o = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinH2O.endf')
return openmc.data.ThermalScattering.from_njoy(
path_h1, path_h2o, temperatures=[293.6, 500.0])
@pytest.fixture(scope='module')
def hzrh():
"""H in ZrH thermal scattering data."""
filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinZrH.endf')
return openmc.data.ThermalScattering.from_endf(filename)
@pytest.fixture(scope='module')
def hzrh_njoy():
"""H in ZrH generated using NJOY."""
path_h1 = os.path.join(_ENDF_DATA, 'neutrons', 'n-001_H_001.endf')
path_hzrh = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinZrH.endf')
with_endf_data = openmc.data.ThermalScattering.from_njoy(
path_h1, path_hzrh, temperatures=[296.0], iwt=0
)
without_endf_data = openmc.data.ThermalScattering.from_njoy(
path_h1, path_hzrh, temperatures=[296.0], use_endf_data=False, iwt=1
)
return with_endf_data, without_endf_data
@pytest.fixture(scope='module')
def sio2():
"""SiO2 thermal scattering data."""
filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-SiO2.endf')
return openmc.data.ThermalScattering.from_endf(filename)
def test_h2o_attributes(h2o):
assert h2o.name == 'c_H_in_H2O'
assert h2o.nuclides == ['H1']
assert h2o.secondary_mode == 'skewed'
assert h2o.temperatures == ['294K']
assert h2o.atomic_weight_ratio == pytest.approx(0.999167)
assert h2o.energy_max == pytest.approx(4.46)
assert isinstance(repr(h2o), str)
def test_h2o_xs(h2o):
assert not h2o.elastic_xs
for temperature, func in h2o.inelastic_xs.items():
assert not h2o.elastic
for temperature, func in h2o.inelastic.xs.items():
assert temperature.endswith('K')
assert isinstance(func, Callable)
@ -53,36 +85,157 @@ def test_h2o_xs(h2o):
def test_graphite_attributes(graphite):
assert graphite.name == 'c_Graphite'
assert graphite.nuclides == ['C0', 'C12', 'C13']
assert graphite.secondary_mode == 'skewed'
assert graphite.temperatures == ['296K']
assert graphite.atomic_weight_ratio == pytest.approx(11.898)
assert graphite.energy_max == pytest.approx(4.46)
def test_graphite_xs(graphite):
for temperature, func in graphite.elastic_xs.items():
for temperature, func in graphite.elastic.xs.items():
assert temperature.endswith('K')
assert isinstance(func, openmc.data.CoherentElastic)
for temperature, func in graphite.inelastic_xs.items():
for temperature, func in graphite.inelastic.xs.items():
assert temperature.endswith('K')
assert isinstance(func, Callable)
elastic = graphite.elastic_xs['296K']
assert elastic([1e-3, 1.0]) == pytest.approx([13.47464936, 0.62590156])
elastic = graphite.elastic.xs['296K']
assert elastic([1e-3, 1.0]) == pytest.approx([0.0, 0.62586153])
def test_export_to_hdf5(tmpdir, h2o_njoy, graphite):
def test_graphite_njoy():
path_c0 = os.path.join(_ENDF_DATA, 'neutrons', 'n-006_C_000.endf')
path_gr = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-graphite.endf')
graphite = openmc.data.ThermalScattering.from_njoy(
path_c0, path_gr, temperatures=[296.0])
assert graphite.nuclides == ['C0', 'C12', 'C13']
assert graphite.atomic_weight_ratio == pytest.approx(11.898)
assert graphite.energy_max == pytest.approx(2.02)
assert graphite.temperatures == ['296K']
def test_export_to_hdf5(tmpdir, h2o_njoy, hzrh_njoy, graphite):
filename = str(tmpdir.join('water.h5'))
h2o_njoy.export_to_hdf5(filename)
assert os.path.exists(filename)
# Graphite covers export of coherent elastic data
filename = str(tmpdir.join('graphite.h5'))
graphite.export_to_hdf5(filename)
assert os.path.exists(filename)
# H in ZrH covers export of incoherent elastic data, and incoherent
# inelastic angle-energy distributions
filename = str(tmpdir.join('hzrh.h5'))
hzrh_njoy[0].export_to_hdf5(filename)
assert os.path.exists(filename)
hzrh_njoy[1].export_to_hdf5(filename, 'w')
assert os.path.exists(filename)
def test_continuous_dist(h2o_njoy):
for temperature, dist in h2o_njoy.inelastic_dist.items():
for temperature, dist in h2o_njoy.inelastic.distribution.items():
assert temperature.endswith('K')
assert isinstance(dist, openmc.data.CorrelatedAngleEnergy)
assert isinstance(dist, openmc.data.IncoherentInelasticAE)
def test_h2o_endf():
filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinH2O.endf')
h2o = openmc.data.ThermalScattering.from_endf(filename)
assert not h2o.elastic
assert h2o.atomic_weight_ratio == pytest.approx(0.99917)
assert h2o.energy_max == pytest.approx(3.99993)
assert h2o.temperatures == ['294K', '350K', '400K', '450K', '500K', '550K',
'600K', '650K', '800K']
def test_hzrh_attributes(hzrh):
assert hzrh.atomic_weight_ratio == pytest.approx(0.99917)
assert hzrh.energy_max == pytest.approx(1.9734)
assert hzrh.temperatures == ['296K', '400K', '500K', '600K', '700K', '800K',
'1000K', '1200K']
def test_hzrh_elastic(hzrh):
rx = hzrh.elastic
for temperature, func in rx.xs.items():
assert temperature.endswith('K')
assert isinstance(func, openmc.data.IncoherentElastic)
xs = rx.xs['296K']
sig_b, W = xs.bound_xs, xs.debye_waller
assert sig_b == pytest.approx(81.98006)
assert W == pytest.approx(8.486993)
for i in range(10):
E = random.uniform(0.0, hzrh.energy_max)
assert xs(E) == pytest.approx(sig_b/2 * ((1 - exp(-4*E*W))/(2*E*W)))
for temperature, dist in rx.distribution.items():
assert temperature.endswith('K')
assert dist.debye_waller > 0.0
def test_hzrh_njoy(hzrh_njoy):
endf, ace = hzrh_njoy
# First check version using ENDF incoherent elastic data
assert endf.atomic_weight_ratio == pytest.approx(0.999167)
assert endf.energy_max == pytest.approx(1.855)
assert endf.temperatures == ['296K']
# Now check version using ACE incoherent elastic data (discretized)
assert ace.atomic_weight_ratio == endf.atomic_weight_ratio
assert ace.energy_max == endf.energy_max
# Cross sections should be about the same (within 1%)
E = np.linspace(1e-5, endf.energy_max)
xs1 = endf.elastic.xs['296K'](E)
xs2 = ace.elastic.xs['296K'](E)
assert xs1 == pytest.approx(xs2, rel=0.01)
# Check discrete incoherent elastic distribution
d = ace.elastic.distribution['296K']
assert np.all((-1.0 <= d.mu_out) & (d.mu_out <= 1.0))
# Check discrete incoherent inelastic distribution
d = endf.inelastic.distribution['296K']
assert d.skewed
assert np.all((-1.0 <= d.mu_out) & (d.mu_out <= 1.0))
assert np.all((0.0 <= d.energy_out) & (d.energy_out < 3*endf.energy_max))
def test_sio2_attributes(sio2):
assert sio2.atomic_weight_ratio == pytest.approx(27.84423)
assert sio2.energy_max == pytest.approx(2.46675)
assert sio2.temperatures == ['294K', '350K', '400K', '500K', '800K',
'1000K', '1200K']
def test_sio2_elastic(sio2):
rx = sio2.elastic
for temperature, func in rx.xs.items():
assert temperature.endswith('K')
assert isinstance(func, openmc.data.CoherentElastic)
xs = rx.xs['294K']
assert len(xs) == 317
assert xs.bragg_edges[0] == pytest.approx(0.000711634)
assert xs.factors[0] == pytest.approx(2.6958e-14)
# Below first bragg edge, cross section should be zero
E = xs.bragg_edges[0] / 2.0
assert xs(E) == 0.0
# Between bragg edges, cross section is P/E where P is the factor
E = (xs.bragg_edges[0] + xs.bragg_edges[1]) / 2.0
P = xs.factors[0]
assert xs(E) == pytest.approx(P / E)
# Check the last Bragg edge
E = 1.1 * xs.bragg_edges[-1]
P = xs.factors[-1]
assert xs(E) == pytest.approx(P / E)
for temperature, dist in rx.distribution.items():
assert temperature.endswith('K')
assert dist.coherent_xs is rx.xs[temperature]
def test_get_thermal_name():
@ -97,5 +250,9 @@ def test_get_thermal_name():
assert f('graphite') == 'c_Graphite'
assert f('D_in_D2O') == 'c_D_in_D2O'
# Not in values, but very close
assert f('hluci') == 'c_H_in_C5O2H8'
assert f('ortho_d') == 'c_ortho_D'
# Names that don't remotely match anything
assert f('boogie_monster') == 'c_boogie_monster'

12
tests/unit_tests/test_deplete_chain.py
View file

@ -231,3 +231,15 @@ def test_getitem():
assert "NucA" == chain["NucA"]
assert "NucB" == chain["NucB"]
assert "NucC" == chain["NucC"]
def test_set_fiss_q():
"""Make sure new fission q values can be set on the chain"""
new_q = {"U235": 2.0E8, "U238": 2.0E8, "U234": 5.0E7}
chain_file = Path(__file__).parents[1] / "chain_simple.xml"
mod_chain = Chain.from_xml(chain_file, new_q)
for name, q in new_q.items():
chain_nuc = mod_chain[name]
for rx in chain_nuc.reactions:
if rx.type == 'fission':
assert rx.Q == q

13
tests/unit_tests/test_deplete_operator.py
View file

@ -68,3 +68,16 @@ def test_operator_init(bare_xs):
'decay_modes', 'yield_data', 'yield_energies',
]:
assert getattr(act_nuc, prop) == getattr(ref_nuc, prop), prop
def test_operator_fiss_q():
"""Make sure fission q values can be set"""
new_q = {"U235": 2.0E8, "U238": 2.0E8, "U234": 5.0E7}
chain_file = Path(__file__).parents[1] / "chain_simple.xml"
operator = BareDepleteOperator(chain_file=chain_file, fission_q=new_q)
mod_chain = operator.chain
for name, q in new_q.items():
chain_nuc = mod_chain[name]
for rx in chain_nuc.reactions:
if rx.type == 'fission':
assert rx.Q == q

2
tools/ci/download-xs.sh
View file

@ -3,7 +3,7 @@ set -ex
# Download HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget -q -O - https://anl.box.com/shared/static/u1g3n8iai0u1n5f6ev3pg2j3ff941bqa.xz | tar -C $HOME -xJ
wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ
fi
# Download ENDF/B-VII.1 distribution