diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f08124811..8d6e622829 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 957bfa3a3c..ef61604cb4 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -215,11 +215,9 @@ Thermal Neutron Scattering Data **//** :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 **//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`. **//inelastic/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 -`. +: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 + `. + +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: -------------------- diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 5cf74bce0a..92cb785481 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -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 -------------- diff --git a/docs/source/usersguide/cross_sections.rst b/docs/source/usersguide/cross_sections.rst index e6a79815f9..636b44db46 100644 --- a/docs/source/usersguide/cross_sections.rst +++ b/docs/source/usersguide/cross_sections.rst @@ -11,7 +11,7 @@ or multi-group mode. In continuous-energy mode, OpenMC uses a native `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 -`_ includes the windowed +`_ includes the windowed multipole library, so if you are using this library, the windowed multipole data will already be available to you. diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 6b2e6f63b7..394f90740f 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -26,7 +26,7 @@ constexpr bool VERSION_DEV {true}; constexpr std::array 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 VERSION_STATEPOINT {17, 0}; @@ -36,7 +36,6 @@ constexpr std::array VERSION_SUMMARY {6, 0}; constexpr std::array VERSION_VOLUME {1, 0}; constexpr std::array VERSION_VOXEL {2, 0}; constexpr std::array VERSION_MGXS_LIBRARY {1, 0}; -constexpr char VERSION_MULTIPOLE[] {"v0.2"}; // ============================================================================ // ADJUSTABLE PARAMETERS diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 7ec078e93f..27958db5c5 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -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 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& 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& bragg_edges() const { return bragg_edges_; } + const std::vector& factors() const { return factors_; } private: std::vector bragg_edges_; //!< Bragg edges in [eV] std::vector 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 diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 8ba18d3d15..603d8c414b 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -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& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index a898db67c8..0217c7e776 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -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 { diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h index 8509b07143..eb90bd8a04 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -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] diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h new file mode 100644 index 0000000000..7b876334f8 --- /dev/null +++ b/include/openmc/secondary_thermal.h @@ -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 +#include "xtensor/xtensor.hpp" + +#include + +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& 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& energy_; //!< Energies at which cosines are tabulated + xt::xtensor 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& 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& energy_; //!< Incident energies + xt::xtensor energy_out_; //!< Outgoing energies for each incident energy + xt::xtensor 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 e_out; //!< Outgoing energies + xt::xtensor e_out_pdf; //!< Probability density function + xt::xtensor e_out_cdf; //!< Cumulative distribution function + xt::xtensor mu; //!< Equiprobable angles at each outgoing energy + }; + + std::vector energy_; //!< Incident energies + std::vector distribution_; //!< Secondary angle-energy at + //!< each incident energy +}; + + +} // namespace openmc + +#endif // OPENMC_SECONDARY_THERMAL_H diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index e895a17bdc..e430c75c45 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -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_; } diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index b737f005b3..a084c18884 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -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 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 e_out; //!< Outgoing energies - xt::xtensor e_out_pdf; //!< Probability density function - xt::xtensor e_out_cdf; //!< Cumulative distribution function - xt::xtensor mu; //!< Equiprobable angles at each outgoing energy + struct Reaction { + // Default constructor + Reaction() { } + + // Data members + std::unique_ptr xs; //!< Cross section + std::unique_ptr 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 inelastic_e_in_; //!< incoming E grid for inelastic - std::vector inelastic_sigma_; //!< inelastic scattering cross section - - // The following are used only for equal/skewed distributions - xt::xtensor inelastic_e_out_; - xt::xtensor 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 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 elastic_e_in_; //!< incoming E grid for elastic - std::vector elastic_P_; //!< elastic scattering cross section - xt::xtensor 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 kTs_; //!< temperatures in [eV] (k*T) std::vector nuclides_; //!< Valid nuclides diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index f52c3be755..277c14ca8d 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -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 * diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index d67cc6b266..8cf9bcf5e8 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -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): diff --git a/openmc/data/correlated.py b/openmc/data/correlated.py index b760f2e2fd..567d23abd4 100644 --- a/openmc/data/correlated.py +++ b/openmc/data/correlated.py @@ -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, diff --git a/openmc/data/decay.py b/openmc/data/decay.py index bbb059f4fa..6a8ad39f09 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -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: diff --git a/openmc/data/function.py b/openmc/data/function.py index c40e2a699c..9175e5e3bf 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -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 ---------- diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index ddc1efb1c4..c8d630851b 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -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. diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index c5abc73c34..05b6b8d2d6 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -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 diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 27ed7aa291..79cad32823 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -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= 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 "".format(self.name) + return "".format(self.name) else: return "" @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 diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py new file mode 100644 index 0000000000..05393e7bb7 --- /dev/null +++ b/openmc/data/thermal_angle_energy.py @@ -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' diff --git a/src/endf.cpp b/src/endf.cpp index 464bc974c2..f01a643113 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -1,6 +1,7 @@ #include "openmc/endf.h" #include // for copy +#include #include // for log, exp #include // for back_inserter #include // for runtime_error @@ -87,6 +88,10 @@ read_function(hid_t group, const char* name) func = std::make_unique(dset); } else if (func_type == "Polynomial") { func = std::make_unique(dset); + } else if (func_type == "CoherentElastic") { + func = std::make_unique(dset); + } else if (func_type == "IncoherentElastic") { + func = std::make_unique(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 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 diff --git a/src/material.cpp b/src/material.cpp index 5667821fb7..ccbbe57a8d 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -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; diff --git a/src/reaction.cpp b/src/reaction.cpp index 222822230d..b733b20283 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -42,8 +42,6 @@ Reaction::Reaction(hid_t group, const std::vector& 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); diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 4f3994ff2f..34680a7735 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -3,6 +3,7 @@ #include // for unique_ptr #include // 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{new Tabulated1D{yield}}; - } else if (temp == "Polynomial") { - yield_ = std::unique_ptr{new Polynomial{yield}}; - } - close_dataset(yield); + yield_ = read_function(group, "yield"); int n; read_attribute(group, "n_distribution", n); diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp new file mode 100644 index 0000000000..696907a17e --- /dev/null +++ b/src/secondary_thermal.cpp @@ -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 // for log, exp + +namespace openmc { + +// Helper function to get index on incident energy grid +void +get_energy_index(const std::vector& 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& 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& 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(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({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 diff --git a/src/thermal.cpp b/src/thermal.cpp index 40bb19e894..5cc395e13c 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -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& 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& 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 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(elastic_.xs.get()); + elastic_.distribution = std::make_unique(*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(elastic_.xs.get()); + if (temp == "incoherent_elastic") { + elastic_.distribution = std::make_unique(dgroup); + } else if (temp == "incoherent_elastic_discrete") { + elastic_.distribution = std::make_unique( + dgroup, xs->x() + ); + } - // Set elastic threshold - threshold_elastic_ = elastic_e_in_.back(); - - // Read angle distribution - if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) { - xt::xarray 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 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(inelastic_.xs.get()); + threshold_inelastic_ = xs->x().back(); - if (secondary_mode != SAB_SECONDARY_CONT) { - // Read energy distribution - xt::xarray 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 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(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({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(dgroup); + } else if (temp == "incoherent_inelastic_discrete") { + inelastic_.distribution = std::make_unique( + 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() diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index 7d7d4e964a..47cbe27b97 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -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' diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 38a5bc4692..07ade9c70f 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -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