diff --git a/CMakeLists.txt b/CMakeLists.txt index 1dee3f90e..4737ff47b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -348,6 +348,7 @@ list(APPEND libopenmc_SOURCES src/message_passing.cpp src/mgxs.cpp src/mgxs_interface.cpp + src/ncrystal_interface.cpp src/nuclide.cpp src/output.cpp src/particle.cpp diff --git a/include/openmc/material.h b/include/openmc/material.h index 5326ad6de..8db9bb628 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -12,13 +12,10 @@ #include "openmc/bremsstrahlung.h" #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr +#include "openmc/ncrystal_interface.h" #include "openmc/particle.h" #include "openmc/vector.h" -#ifdef NCRYSTAL -#include "NCrystal/NCrystal.hh" -#endif - namespace openmc { //============================================================================== @@ -155,25 +152,17 @@ public: //! \return Temperature in [K] double temperature() const; -#ifdef NCRYSTAL //! Get pointer to NCrystal material object //! \return Pointer to NCrystal material object - std::shared_ptr ncrystal_mat() const - { - return ncrystal_mat_; - }; -#endif + const NCrystalMat& ncrystal_mat() const { return ncrystal_mat_; }; //---------------------------------------------------------------------------- // Data - int32_t id_ {C_NONE}; //!< Unique ID - std::string name_; //!< Name of material - vector nuclide_; //!< Indices in nuclides vector - vector element_; //!< Indices in elements vector -#ifdef NCRYSTAL - std::string ncrystal_cfg_; //!< NCrystal configuration string - std::shared_ptr ncrystal_mat_; -#endif + int32_t id_ {C_NONE}; //!< Unique ID + std::string name_; //!< Name of material + vector nuclide_; //!< Indices in nuclides vector + vector element_; //!< Indices in elements vector + NCrystalMat ncrystal_mat_; //!< NCrystal material object xt::xtensor atom_density_; //!< Nuclide atom density in [atom/b-cm] double density_; //!< Total atom density in [atom/b-cm] double density_gpcc_; //!< Total atom density in [g/cm^3] diff --git a/include/openmc/ncrystal_interface.h b/include/openmc/ncrystal_interface.h new file mode 100644 index 000000000..36038d50c --- /dev/null +++ b/include/openmc/ncrystal_interface.h @@ -0,0 +1,91 @@ +#ifndef OPENMC_NCRYSTAL_INTERFACE_H +#define OPENMC_NCRYSTAL_INTERFACE_H + +#ifdef NCRYSTAL +#include "NCrystal/NCRNG.hh" +#include "NCrystal/NCrystal.hh" +#endif + +#include "openmc/particle.h" + +#include // for uint64_t +#include // for numeric_limits +#include + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +extern "C" const bool NCRYSTAL_ENABLED; + +//============================================================================== +// Wrapper class an NCrystal material +//============================================================================== + +class NCrystalMat { +public: + //---------------------------------------------------------------------------- + // Constructors + NCrystalMat() = default; + explicit NCrystalMat(const std::string& cfg); + + //---------------------------------------------------------------------------- + // Methods + +#ifdef NCRYSTAL + //! Return configuration string + std::string cfg() const; + + //! Get cross section from NCrystal material + // + //! \param[in] p Particle object + //! \return Cross section in [b] + double xs(const Particle& p) const; + + // Process scattering event + // + //! \param[in] p Particle object + void scatter(Particle& p) const; + + //! Whether the object holds a valid NCrystal material + operator bool() const; +#else + + //---------------------------------------------------------------------------- + // Trivial methods when compiling without NCRYSTAL + std::string cfg() const + { + return ""; + } + double xs(const Particle& p) const + { + return -1.0; + } + void scatter(Particle& p) const {} + operator bool() const + { + return false; + } +#endif + +private: + //---------------------------------------------------------------------------- + // Data members (only present when compiling with NCrystal support) +#ifdef NCRYSTAL + std::string cfg_; //!< NCrystal configuration string + std::shared_ptr + ptr_; //!< Pointer to NCrystal material object +#endif +}; + +//============================================================================== +// Functions +//============================================================================== + +void ncrystal_update_micro(double xs, NuclideMicroXS& micro); + +} // namespace openmc + +#endif // OPENMC_NCRYSTAL_INTERFACE_H diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 806297c4e..382f0f0dd 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -8,10 +8,6 @@ #include "openmc/reaction.h" #include "openmc/vector.h" -#ifdef NCRYSTAL -#include "NCrystal/NCRNG.hh" -#endif - namespace openmc { //============================================================================== @@ -100,31 +96,11 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p); void sample_secondary_photons(Particle& p, int i_nuclide); -//！Split or Roulette particles based their weight and the lower weight window -// bound. +// ！Split or Roulette particles based their weight and the lower weight window +// bound. //! \param[in] p, particle to be split or rouletted with the weight window. void split_particle(Particle& p); -#ifdef NCRYSTAL -//============================================================================== -// NCrystal wrapper class for the OpenMC random number generator -//============================================================================== - -class NCrystalRNGWrapper : public NCrystal::RNGStream { -public: - constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} - -protected: - double actualGenerate() override - { - return std::max( - std::numeric_limits::min(), prn(openmc_seed_)); - } -private: - uint64_t* openmc_seed_; -}; -#endif - } // namespace openmc #endif // OPENMC_PHYSICS_H diff --git a/include/openmc/settings.h b/include/openmc/settings.h index fb5f8387a..9cd430deb 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -20,8 +20,6 @@ namespace openmc { // Global variable declarations //============================================================================== -extern "C" const bool NCRYSTAL_ENABLED; - namespace settings { // Boolean flags @@ -93,6 +91,8 @@ extern int n_log_bins; //!< number of bins for logarithmic energy grid extern int n_batches; //!< number of (inactive+active) batches extern int n_max_batches; //!< Maximum number of batches extern int max_tracks; //!< Maximum number of particle tracks written to file +extern double + ncrystal_max_energy; //!< Energy in eV to switch between NCrystal and ENDF extern ResScatMethod res_scat_method; //!< resonance upscattering method extern double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering extern double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering @@ -124,10 +124,6 @@ extern int trigger_batch_interval; //!< Batch interval for triggers extern "C" int verbosity; //!< How verbose to make output extern double weight_cutoff; //!< Weight cutoff for Russian roulette extern double weight_survive; //!< Survival weight after Russian roulette -#ifdef NCRYSTAL -extern double - ncrystal_max_energy; //!< Energy in eV to switch between NCrystal and ENDF -#endif } // namespace settings //============================================================================== diff --git a/src/material.cpp b/src/material.cpp index 25014c448..2d55d6a04 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -60,13 +60,12 @@ Material::Material(pugi::xml_node node) name_ = get_node_value(node, "name"); } -#ifdef NCRYSTAL if (check_for_node(node, "cfg")) { - ncrystal_cfg_ = get_node_value(node, "cfg"); - write_message(5, "NCrystal config string for material #{}: '{}'", this->id(), ncrystal_cfg_); - ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); + auto cfg = get_node_value(node, "cfg"); + write_message( + 5, "NCrystal config string for material #{}: '{}'", this->id(), cfg); + ncrystal_mat_ = NCrystalMat(cfg); } -#endif if (check_for_node(node, "depletable")) { depletable_ = get_node_value_bool(node, "depletable"); @@ -800,19 +799,11 @@ void Material::calculate_neutron_xs(Particle& p) const // Initialize position in i_sab_nuclides int j = 0; -#ifdef NCRYSTAL - double ncrystal_xs = -1; - + // Calculate NCrystal cross section + double ncrystal_xs = -1.0; if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy) { - // Calculate scattering XS per atom with NCrystal, only once per material - NCrystal::CachePtr dummy_cache; - auto nc_energy = NCrystal::NeutronEnergy {p.E()}; - ncrystal_xs = - ncrystal_mat_ - ->crossSection(dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z}) - .get(); + ncrystal_xs = ncrystal_mat_.xs(p); } -#endif // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { @@ -852,26 +843,16 @@ void Material::calculate_neutron_xs(Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide -#ifdef NCRYSTAL auto& micro {p.neutron_xs(i_nuclide)}; -#else - const auto& micro {p.neutron_xs(i_nuclide)}; -#endif if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT || i_sab != micro.index_sab || sab_frac != micro.sab_frac) { data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p); -#ifdef NCRYSTAL + + // If NCrystal is being used, update micro cross section cache if (ncrystal_xs >= 0.0) { - if (micro.thermal > 0 || micro.thermal_elastic > 0) { - fatal_error("S(a,b) treatment and NCrystal are not compatible."); - } data::nuclides[i_nuclide]->calculate_elastic_xs(p); - // remove free atom cross section - // and replace it by scattering cross section per atom from NCrystal - micro.total = micro.total - micro.elastic + ncrystal_xs; - micro.elastic = ncrystal_xs; + ncrystal_update_micro(ncrystal_xs, micro); } -#endif } // ====================================================================== diff --git a/src/ncrystal_interface.cpp b/src/ncrystal_interface.cpp new file mode 100644 index 000000000..b39f62d90 --- /dev/null +++ b/src/ncrystal_interface.cpp @@ -0,0 +1,108 @@ +#include "openmc/ncrystal_interface.h" + +#include "openmc/error.h" +#include "openmc/material.h" +#include "openmc/random_lcg.h" + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +#ifdef NCRYSTAL +const bool NCRYSTAL_ENABLED = true; +#else +const bool NCRYSTAL_ENABLED = false; +#endif + +//============================================================================== +// NCrystal wrapper class for the OpenMC random number generator +//============================================================================== + +#ifdef NCRYSTAL +class NCrystalRNGWrapper : public NCrystal::RNGStream { +public: + constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} + +protected: + double actualGenerate() override + { + return std::max( + std::numeric_limits::min(), prn(openmc_seed_)); + } + +private: + uint64_t* openmc_seed_; +}; +#endif + +//============================================================================== +// NCrystal implementation +//============================================================================== + +NCrystalMat::NCrystalMat(const std::string& cfg) +{ +#ifdef NCRYSTAL + cfg_ = cfg; + ptr_ = NCrystal::FactImpl::createScatter(cfg); +#else + fatal_error("Your build of OpenMC does not support NCrystal materials."); +#endif +} + +#ifdef NCRYSTAL +std::string NCrystalMat::cfg() const +{ + return cfg_; +} + +double NCrystalMat::xs(const Particle& p) const +{ + // Calculate scattering XS per atom with NCrystal, only once per material + NCrystal::CachePtr dummy_cache; + auto nc_energy = NCrystal::NeutronEnergy {p.E()}; + return ptr_->crossSection(dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z}) + .get(); +} + +void NCrystalMat::scatter(Particle& p) const +{ + NCrystalRNGWrapper rng(p.current_seed()); // Initialize RNG + // create a cache pointer for multi thread physics + NCrystal::CachePtr dummy_cache; + auto nc_energy = NCrystal::NeutronEnergy {p.E()}; + auto outcome = ptr_->sampleScatter( + dummy_cache, rng, nc_energy, {p.u().x, p.u().y, p.u().z}); + + // Modify attributes of particle + p.E() = outcome.ekin.get(); + Direction u_old {p.u()}; + p.u() = + Direction(outcome.direction[0], outcome.direction[1], outcome.direction[2]); + p.mu() = u_old.dot(p.u()); + p.event_mt() = ELASTIC; +} + +NCrystalMat::operator bool() const +{ + return ptr_.get(); +} +#endif + +//============================================================================== +// Functions +//============================================================================== + +void ncrystal_update_micro(double xs, NuclideMicroXS& micro) +{ + if (micro.thermal > 0 || micro.thermal_elastic > 0) { + fatal_error("S(a,b) treatment and NCrystal are not compatible."); + } + // remove free atom cross section + // and replace it by scattering cross section per atom from NCrystal + micro.total = micro.total - micro.elastic + xs; + micro.elastic = xs; +} + +} // namespace openmc diff --git a/src/physics.cpp b/src/physics.cpp index c1ef282f7..46bf37aa9 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -10,6 +10,7 @@ #include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" +#include "openmc/ncrystal_interface.h" #include "openmc/nuclide.h" #include "openmc/photon.h" #include "openmc/physics_common.h" @@ -138,32 +139,15 @@ void sample_neutron_reaction(Particle& p) if (!p.alive()) return; - // Sample a scattering reaction and determine the secondary energy of the - // exiting neutron - -#ifdef NCRYSTAL - if (model::materials[p.material()]->ncrystal_mat() && - p.E() < settings::ncrystal_max_energy) { - NCrystalRNGWrapper rng(p.current_seed()); // Initialize RNG - // create a cache pointer for multi thread physics - NCrystal::CachePtr dummy_cache; - auto nc_energy = NCrystal::NeutronEnergy {p.E()}; - auto outcome = - model::materials[p.material()]->ncrystal_mat()->sampleScatter( - dummy_cache, rng, nc_energy, {p.u().x, p.u().y, p.u().z}); - p.E_last() = p.E(); - p.E() = outcome.ekin.get(); - Direction u_old {p.u()}; - p.u() = Direction( - outcome.direction[0], outcome.direction[1], outcome.direction[2]); - p.mu() = u_old.dot(p.u()); - p.event_mt() = ELASTIC; + // Sample a scattering reaction and determine the secondary energy of the + // exiting neutron + const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); + if (ncrystal_mat && p.E() < settings::ncrystal_max_energy) { + ncrystal_mat.scatter(p); } else { scatter(p, i_nuclide); } -#else - scatter(p, i_nuclide); -#endif + // Advance URR seed stream 'N' times after energy changes if (p.E() != p.E_last()) { p.stream() = STREAM_URR_PTABLE; diff --git a/src/settings.cpp b/src/settings.cpp index 6f6e4e090..9174d4114 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -37,13 +37,6 @@ namespace openmc { // Global variables //============================================================================== -#ifdef NCRYSTAL -const bool NCRYSTAL_ENABLED = true; -#else -const bool NCRYSTAL_ENABLED = false; -#endif - - namespace settings { // Default values for boolean flags @@ -106,6 +99,7 @@ int n_batches; int n_max_batches; int max_splits {1000}; int max_tracks {1000}; +double ncrystal_max_energy {5.0}; ResScatMethod res_scat_method {ResScatMethod::rvs}; double res_scat_energy_min {0.01}; double res_scat_energy_max {1000.0}; @@ -128,10 +122,6 @@ int verbosity {7}; double weight_cutoff {0.25}; double weight_survive {1.0}; -#ifdef NCRYSTAL -double ncrystal_max_energy {5.0}; -#endif - } // namespace settings //==============================================================================