From 0e0cd0a22d059159a8799906f57d12cec99c1a89 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Mon, 13 Sep 2021 01:25:36 -0300 Subject: [PATCH] Added hooks for NCrystal --- CMakeLists.txt | 16 +++++++++ include/openmc/material.h | 8 +++++ include/openmc/particle_data.h | 3 ++ include/openmc/physics.h | 27 ++++++++++++++++ include/openmc/settings.h | 3 ++ openmc/material.py | 59 ++++++++++++++++++++++++++++++++++ src/material.cpp | 36 ++++++++++++++++++++- src/physics.cpp | 21 +++++++++++- src/settings.cpp | 4 +++ 9 files changed, 175 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 521b46cfc6..7be54dc3e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,7 @@ option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags" option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" OFF) option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF) option(OPENMC_USE_MPI "Enable MPI" OFF) +option(OPENMC_USE_NCRYSTAL "Enable support for NCrystal scattering" OFF) # Warnings for deprecated options foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh") @@ -98,6 +99,14 @@ macro(find_package_write_status pkg) endif() endmacro() +#=============================================================================== +# NCrystal Scattering Support +#=============================================================================== + +if(OPENMC_USE_NCRYSTAL) + find_package(NCrystal REQUIRED PATH_SUFFIXES ./)#fixme +endif() + #=============================================================================== # DAGMC Geometry Support - need DAGMC/MOAB #=============================================================================== @@ -482,6 +491,13 @@ if (OPENMC_USE_MPI) target_link_libraries(libopenmc MPI::MPI_CXX) endif() +if(OPENMC_USE_NCRYSTAL) + target_compile_definitions(libopenmc PRIVATE NCRYSTAL) + target_link_libraries(libopenmc NCrystal::NCrystal) +endif() + + + #=============================================================================== # openmc executable #=============================================================================== diff --git a/include/openmc/material.h b/include/openmc/material.h index b251a3ca85..55f8ae20d1 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -15,6 +15,10 @@ #include "openmc/particle.h" #include "openmc/vector.h" +#ifdef NCRYSTAL +#include "NCrystal/NCrystal.hh" +#endif + namespace openmc { //============================================================================== @@ -157,6 +161,10 @@ public: std::string name_; //!< Name of material vector nuclide_; //!< Indices in nuclides vector vector element_; //!< Indices in elements vector +#ifdef NCRYSTAL + std::string cfg_; //!< NCrystal configuration string + std::shared_ptr m_NCrystal_mat_; +#endif 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/particle_data.h b/include/openmc/particle_data.h index d3d00571a8..06ca5990a5 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -162,6 +162,9 @@ struct MacroXS { double fission; //!< macroscopic fission xs double nu_fission; //!< macroscopic production xs double photon_prod; //!< macroscopic photon production xs +#ifdef NCRYSTAL + double NCrystal_XS; //!< macroscopic cross section of processes handled by NCrystal +#endif // Photon cross sections double coherent; //!< macroscopic coherent xs diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 262b3a8841..a1db072914 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -8,6 +8,10 @@ #include "openmc/reaction.h" #include "openmc/vector.h" +#ifdef NCRYSTAL +#include "NCrystal/NCRNG.hh" +#endif + namespace openmc { //============================================================================== @@ -101,6 +105,29 @@ void sample_secondary_photons(Particle& p, int i_nuclide); //! \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 NcrystalRNG_Wrapper : public NCrystal::RNGStream { + uint64_t* openmc_seed_; +public: + constexpr NcrystalRNG_Wrapper(uint64_t* s) noexcept : openmc_seed_(s) {} + //Can be cheaply created on the stack just before being used in calls to + //ProcImpl::Scatter objects, like: + // + // RNG_Wrapper rng(seed); + // +protected: + //double actualGenerate() override { return prn(openmc_seed_); } + double actualGenerate() override { + return std::max( std::numeric_limits::min(), prn(openmc_seed_) ); + } + +}; +#endif + } // namespace openmc #endif // OPENMC_PHYSICS_H diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 1e061b235b..ba1a30e307 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -120,6 +120,9 @@ 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/openmc/material.py b/openmc/material.py index e99eec52c4..1fe62961d9 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -118,6 +118,7 @@ class Material(IDManagerMixin): self._volume = None self._atoms = {} self._isotropic = [] + self._NCrystal_cfg = None # A list of tuples (nuclide, percent, percent type) self._nuclides = [] @@ -140,6 +141,8 @@ class Material(IDManagerMixin): string += '{: <16}\n'.format('\tS(a,b) Tables') + string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._NCrystal_cfg) + for sab in self._sab: string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab) @@ -331,6 +334,48 @@ class Material(IDManagerMixin): return material + @classmethod + def from_NCrystal(cls, cfg): + """Create material from NCrystal configuration string + Density is set from the NCrystal value, + and material temperature from the configuration string. + + Parameters + ---------- + cfg : str + NCrystal configuration string + + Returns + ------- + openmc.Material + Material instance + + """ + + try: + import NCrystal + except ImportError: + raise SystemExit("ERROR: NCrystal Python module is required.") + + NC_mat = NCrystal.createInfo(cfg) + NC_comp = NC_mat.getComposition() + + # Create the Material + material = cls() + + for frac, atom in NC_comp: + if (atom.A() == 0): + material.add_element(atom.displayLabel(), frac, 'ao') + else: + material.add_nuclide(atom.displayLabel(), frac, 'ao') + + material._NCrystal_cfg = cfg + material._density_units = "g/cm3" + material._density = NC_mat.getDensity() + material.temperature = NC_mat.getTemperature() + + return material + def add_volume_information(self, volume_calc): """Add volume information to a material. @@ -405,6 +450,9 @@ class Material(IDManagerMixin): 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) + if self._NCrystal_cfg is not None: + raise ValueError("Cannot add nuclides to NCrystal material") + # If nuclide name doesn't look valid, give a warning try: Z, _, _ = openmc.data.zam(nuclide) @@ -609,6 +657,9 @@ class Material(IDManagerMixin): raise ValueError("Element name should be given by the " "element's symbol or name, e.g., 'Zr', 'zirconium'") + if self._NCrystal_cfg is not None: + raise ValueError("Cannot add elements to NCrystal material") + # Allow for element identifier to be given as a symbol or name if len(element) > 2: el = element.lower() @@ -1135,6 +1186,14 @@ class Material(IDManagerMixin): if self._volume: element.set("volume", str(self._volume)) + if self._NCrystal_cfg: + if self._sab: + raise ValueError("NCrystal materials are not compatible with S(a,b).") + if self._macroscopic is not None: + raise ValueError("NCrystal materials are not compatible macroscopic cross sections.") + + element.set("cfg", str(self._NCrystal_cfg)) + # Create temperature XML subelement if self.temperature is not None: element.set("temperature", str(self.temperature)) diff --git a/src/material.cpp b/src/material.cpp index 30dfa5ed50..d7fae83617 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -60,6 +60,14 @@ Material::Material(pugi::xml_node node) name_ = get_node_value(node, "name"); } +#ifdef NCRYSTAL + if (check_for_node(node, "cfg")) { + cfg_ = get_node_value(node, "cfg"); + write_message(5, "NCrystal config string: >>{}<< ", cfg_); + m_NCrystal_mat_ = NCrystal::FactImpl::createScatter(cfg_); + } +#endif + if (check_for_node(node, "depletable")) { depletable_ = get_node_value_bool(node, "depletable"); } @@ -792,6 +800,17 @@ void Material::calculate_neutron_xs(Particle& p) const // Initialize position in i_sab_nuclides int j = 0; +#ifdef NCRYSTAL + double ncrystal_xs = -1; + + if (m_NCrystal_mat_ != nullptr && p.E() < settings::ncrystal_max_energy){ + // Calculate scattering XS per atom with NCrystal, only once per material + NCrystal::CachePtr dummyCache; + auto nc_energy = NCrystal::NeutronEnergy{p.E()}; + ncrystal_xs = m_NCrystal_mat_->crossSection( dummyCache, nc_energy, {p.u().x, p.u().y, p.u().z}).get(); + } +#endif + // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { // ====================================================================== @@ -830,10 +849,25 @@ void Material::calculate_neutron_xs(Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {p.neutron_xs(i_nuclide)}; +#ifndef NCRYSTAL + const +#endif + auto& micro {p.neutron_xs(i_nuclide)}; 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_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; + } +#endif } // ====================================================================== diff --git a/src/physics.cpp b/src/physics.cpp index bcdcf7ecf9..9eff4aed7e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -140,8 +140,27 @@ void sample_neutron_reaction(Particle& p) // Sample a scattering reaction and determine the secondary energy of the // exiting neutron - scatter(p, i_nuclide); +#ifdef NCRYSTAL + if (model::materials[p.material_]->m_NCrystal_mat_ != nullptr && p.E() < settings::ncrystal_max_energy){ + + NcrystalRNG_Wrapper rng(p.current_seed()); // Initialize RNG + //create a cache pointer for multi thread physics + NCrystal::CachePtr dummyCache;//fixme: avoid recreating here (triggers malloc) + auto nc_energy = NCrystal::NeutronEnergy{p.E()}; + auto outcome = model::materials[p.material_]->m_NCrystal_mat_->sampleScatter( dummyCache, 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;//fixme: define a new label for NCrystal? + } 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 eb0d4936c7..3f3e970541 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -118,6 +118,10 @@ int verbosity {7}; double weight_cutoff {0.25}; double weight_survive {1.0}; +#ifdef NCRYSTAL +double ncrystal_max_energy{5.0}; +#endif + } // namespace settings //==============================================================================