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 01/42] 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 521b46cfc..7be54dc3e 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 b251a3ca8..55f8ae20d 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 d3d00571a..06ca5990a 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 262b3a884..a1db07291 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 1e061b235..ba1a30e30 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 e99eec52c..1fe62961d 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 30dfa5ed5..d7fae8361 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 bcdcf7ecf..9eff4aed7 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 eb0d4936c..3f3e97054 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 //============================================================================== From 0f33e0b11c6add029f354eb3e5c89b002eec7500 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Tue, 13 Sep 2022 14:22:04 +0200 Subject: [PATCH 02/42] Add temporary requirement of natural elements in NCrystal materials --- openmc/material.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 1fe62961d..0a293eb0f 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -364,10 +364,9 @@ class Material(IDManagerMixin): 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') + if not atom.isNaturalElement(): + raise ValueError('NCrystal-OpenMC interface only works with natural elements for now.') + material.add_element(atom.elementName(), frac, 'ao') material._NCrystal_cfg = cfg material._density_units = "g/cm3" From 9c54ac0b6b53336f834b5f868aa85e3e140239f8 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Tue, 13 Sep 2022 16:07:50 +0200 Subject: [PATCH 03/42] Fix access to private attributes --- include/openmc/material.h | 6 ++++++ src/physics.cpp | 9 ++++----- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index 55f8ae20d..eef557088 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -155,6 +155,12 @@ public: //! \return Temperature in [K] double temperature() const; +#ifdef NCRYSTAL + //! Gwet pointer to NCrystal material object + //! \return Pointer to NCrystal material object + std::shared_ptr m_NCrystal_mat() const {return m_NCrystal_mat_}; +#endif + //---------------------------------------------------------------------------- // Data int32_t id_ {C_NONE}; //!< Unique ID diff --git a/src/physics.cpp b/src/physics.cpp index 9eff4aed7..948307521 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -142,19 +142,18 @@ void sample_neutron_reaction(Particle& p) // exiting neutron #ifdef NCRYSTAL - if (model::materials[p.material_]->m_NCrystal_mat_ != nullptr && p.E() < settings::ncrystal_max_energy){ - + 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}); + 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? + p.mu() = u_old.dot(p.u()); + p.event_mt() = ELASTIC;//fixme: define a new label for NCrystal? } else { scatter(p, i_nuclide); } From 70d5c76ebd70cf43349423f17028c3cb1dfe1daa Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Tue, 13 Sep 2022 16:10:05 +0200 Subject: [PATCH 04/42] Fixed material.h --- include/openmc/material.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index eef557088..10e64773e 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -158,7 +158,7 @@ public: #ifdef NCRYSTAL //! Gwet pointer to NCrystal material object //! \return Pointer to NCrystal material object - std::shared_ptr m_NCrystal_mat() const {return m_NCrystal_mat_}; + std::shared_ptr m_NCrystal_mat() const {return m_NCrystal_mat_; }; #endif //---------------------------------------------------------------------------- From 29a10eafd3db601ffa9a7e7bf882431c150ee1b2 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Fri, 16 Sep 2022 09:03:09 -0300 Subject: [PATCH 05/42] Removed MacroNCrystalXS This was left from a previous implementation. Now it is not needed because NCrystal just takes over all neutron scattering events below a given energy threshold. --- include/openmc/particle_data.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 06ca5990a..d3d00571a 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -162,9 +162,6 @@ 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 From 2e19260dea9176f03f2d5d80be19c8c819d339d3 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Mon, 26 Sep 2022 11:19:54 +0200 Subject: [PATCH 06/42] Implement recommended changes --- CMakeLists.txt | 2 +- include/openmc/material.h | 8 ++++---- include/openmc/physics.h | 10 ++-------- openmc/material.py | 31 ++++++++++++++----------------- src/material.cpp | 16 +++++++++------- src/physics.cpp | 10 +++++----- 6 files changed, 35 insertions(+), 42 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7be54dc3e..7080accfd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -104,7 +104,7 @@ endmacro() #=============================================================================== if(OPENMC_USE_NCRYSTAL) - find_package(NCrystal REQUIRED PATH_SUFFIXES ./)#fixme + find_package(NCrystal REQUIRED) endif() #=============================================================================== diff --git a/include/openmc/material.h b/include/openmc/material.h index 10e64773e..30d8ffb7d 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -156,9 +156,9 @@ public: double temperature() const; #ifdef NCRYSTAL - //! Gwet pointer to NCrystal material object + //! Get pointer to NCrystal material object //! \return Pointer to NCrystal material object - std::shared_ptr m_NCrystal_mat() const {return m_NCrystal_mat_; }; + std::shared_ptr ncrystal_mat() const {return ncrystal_mat_; }; #endif //---------------------------------------------------------------------------- @@ -168,8 +168,8 @@ public: 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_; + std::string ncrystal_cfg_; //!< NCrystal configuration string + std::shared_ptr ncrystal_mat_; #endif xt::xtensor atom_density_; //!< Nuclide atom density in [atom/b-cm] double density_; //!< Total atom density in [atom/b-cm] diff --git a/include/openmc/physics.h b/include/openmc/physics.h index a1db07291..109f7caff 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -110,17 +110,11 @@ void split_particle(Particle& p); // NCrystal wrapper class for the OpenMC random number generator //============================================================================== -class NcrystalRNG_Wrapper : public NCrystal::RNGStream { +class NCrystalRNGWrapper : 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); - // + constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} protected: - //double actualGenerate() override { return prn(openmc_seed_); } double actualGenerate() override { return std::max( std::numeric_limits::min(), prn(openmc_seed_) ); } diff --git a/openmc/material.py b/openmc/material.py index 0a293eb0f..15d1f8ac9 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -118,7 +118,7 @@ class Material(IDManagerMixin): self._volume = None self._atoms = {} self._isotropic = [] - self._NCrystal_cfg = None + self._ncrystal_cfg = None # A list of tuples (nuclide, percent, percent type) self._nuclides = [] @@ -141,7 +141,7 @@ class Material(IDManagerMixin): string += '{: <16}\n'.format('\tS(a,b) Tables') - string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._NCrystal_cfg) + string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg) for sab in self._sab: string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab) @@ -335,7 +335,7 @@ class Material(IDManagerMixin): return material @classmethod - def from_NCrystal(cls, cfg): + 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. @@ -352,26 +352,23 @@ class Material(IDManagerMixin): """ - try: - import NCrystal - except ImportError: - raise SystemExit("ERROR: NCrystal Python module is required.") + import NCrystal - NC_mat = NCrystal.createInfo(cfg) - NC_comp = NC_mat.getComposition() + nc_mat = NCrystal.createInfo(cfg) + nc_comp = nc_mat.getComposition() # Create the Material material = cls() - for frac, atom in NC_comp: + for frac, atom in nc_comp: if not atom.isNaturalElement(): raise ValueError('NCrystal-OpenMC interface only works with natural elements for now.') material.add_element(atom.elementName(), frac, 'ao') - material._NCrystal_cfg = cfg + material._ncrystal_cfg = cfg material._density_units = "g/cm3" - material._density = NC_mat.getDensity() - material.temperature = NC_mat.getTemperature() + material._density = nc_mat.getDensity() + material.temperature = nc_mat.getTemperature() return material @@ -449,7 +446,7 @@ class Material(IDManagerMixin): 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) - if self._NCrystal_cfg is not None: + 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 @@ -656,7 +653,7 @@ 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: + 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 @@ -1185,13 +1182,13 @@ class Material(IDManagerMixin): if self._volume: element.set("volume", str(self._volume)) - if self._NCrystal_cfg: + 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)) + element.set("cfg", str(self._ncrystal_cfg)) # Create temperature XML subelement if self.temperature is not None: diff --git a/src/material.cpp b/src/material.cpp index d7fae8361..13acbba1f 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -63,8 +63,8 @@ Material::Material(pugi::xml_node node) #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_); + write_message(5, "NCrystal config string: '{}'", ncrystal_cfg_); + ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); } #endif @@ -803,11 +803,11 @@ void Material::calculate_neutron_xs(Particle& p) const #ifdef NCRYSTAL double ncrystal_xs = -1; - if (m_NCrystal_mat_ != nullptr && p.E() < settings::ncrystal_max_energy){ + if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy){ // Calculate scattering XS per atom with NCrystal, only once per material - NCrystal::CachePtr dummyCache; + NCrystal::CachePtr dummy_cache; 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(); + ncrystal_xs = ncrystal_mat_->crossSection( dummy_cache, nc_energy, {p.u().x, p.u().y, p.u().z}).get(); } #endif @@ -849,8 +849,10 @@ void Material::calculate_neutron_xs(Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide -#ifndef NCRYSTAL - const +#ifdef NCRYSTAL + auto& micro {p.neutron_xs(i_nuclide)}; +#else + const auto& micro {p.neutron_xs(i_nuclide)}; #endif auto& micro {p.neutron_xs(i_nuclide)}; if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT || diff --git a/src/physics.cpp b/src/physics.cpp index 948307521..bad499772 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -142,18 +142,18 @@ void sample_neutron_reaction(Particle& p) // exiting neutron #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 + 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 dummyCache;//fixme: avoid recreating here (triggers malloc) + NCrystal::CachePtr dummy_cache; 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}); + 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;//fixme: define a new label for NCrystal? + p.event_mt() = ELASTIC; } else { scatter(p, i_nuclide); } From 5afcc9e49b60ec1c5e2290d4725bac7a96e1692e Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Mon, 26 Sep 2022 11:30:39 +0200 Subject: [PATCH 07/42] Apply clang-format --- include/openmc/material.h | 13 ++++++++----- include/openmc/physics.h | 17 ++++++++++------- include/openmc/settings.h | 3 ++- src/material.cpp | 33 ++++++++++++++++++--------------- src/physics.cpp | 18 +++++++++++------- src/settings.cpp | 2 +- 6 files changed, 50 insertions(+), 36 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index 30d8ffb7d..47e3e6a16 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -158,15 +158,18 @@ public: #ifdef NCRYSTAL //! Get pointer to NCrystal material object //! \return Pointer to NCrystal material object - std::shared_ptr ncrystal_mat() const {return ncrystal_mat_; }; + std::shared_ptr ncrystal_mat() const + { + return ncrystal_mat_; + }; #endif //---------------------------------------------------------------------------- // 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 + 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_; diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 109f7caff..23995e7af 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -111,14 +111,17 @@ void split_particle(Particle& p); //============================================================================== class NCrystalRNGWrapper : public NCrystal::RNGStream { - uint64_t* openmc_seed_; -public: - constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} -protected: - double actualGenerate() override { - return std::max( std::numeric_limits::min(), prn(openmc_seed_) ); - } + uint64_t* openmc_seed_; +public: + constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} + +protected: + double actualGenerate() override + { + return std::max( + std::numeric_limits::min(), prn(openmc_seed_)); + } }; #endif diff --git a/include/openmc/settings.h b/include/openmc/settings.h index ba1a30e30..9bb447005 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -121,7 +121,8 @@ 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 +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 13acbba1f..e894084e2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -64,7 +64,7 @@ Material::Material(pugi::xml_node node) if (check_for_node(node, "cfg")) { cfg_ = get_node_value(node, "cfg"); write_message(5, "NCrystal config string: '{}'", ncrystal_cfg_); - ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); + ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); } #endif @@ -375,8 +375,8 @@ void Material::finalize() this->init_thermal(); } -// Normalize density -this->normalize_density(); + // Normalize density + this->normalize_density(); } void Material::normalize_density() @@ -803,11 +803,14 @@ void Material::calculate_neutron_xs(Particle& p) const #ifdef NCRYSTAL double ncrystal_xs = -1; - if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy){ + 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(); + 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(); } #endif @@ -859,15 +862,15 @@ void Material::calculate_neutron_xs(Particle& p) const 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; + 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 bad499772..e19640314 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -138,20 +138,24 @@ void sample_neutron_reaction(Particle& p) if (!p.alive()) return; - // Sample a scattering reaction and determine the secondary energy of the - // exiting neutron + // 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){ + 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 + // 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}); + 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.u() = Direction( + outcome.direction[0], outcome.direction[1], outcome.direction[2]); p.mu() = u_old.dot(p.u()); p.event_mt() = ELASTIC; } else { diff --git a/src/settings.cpp b/src/settings.cpp index 3f3e97054..87c24eedd 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -119,7 +119,7 @@ double weight_cutoff {0.25}; double weight_survive {1.0}; #ifdef NCRYSTAL -double ncrystal_max_energy{5.0}; +double ncrystal_max_energy {5.0}; #endif } // namespace settings From fdf763e801b1439da96009fd5dadeadeb93e36f7 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Mon, 26 Sep 2022 12:35:12 +0200 Subject: [PATCH 08/42] Couple of fixes --- src/material.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/material.cpp b/src/material.cpp index e894084e2..9531e9f43 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -62,7 +62,7 @@ Material::Material(pugi::xml_node node) #ifdef NCRYSTAL if (check_for_node(node, "cfg")) { - cfg_ = get_node_value(node, "cfg"); + ncrystal_cfg_ = get_node_value(node, "cfg"); write_message(5, "NCrystal config string: '{}'", ncrystal_cfg_); ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); } @@ -857,7 +857,6 @@ void Material::calculate_neutron_xs(Particle& p) const #else const auto& micro {p.neutron_xs(i_nuclide)}; #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); From 52ee2c37ced8fce859d2ec6a06c6909dd0d04b44 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Thu, 20 Oct 2022 11:36:22 +0200 Subject: [PATCH 09/42] Print where NCrystal was found as suggested by @paulromano Co-authored-by: Paul Romano --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7080accfd..2815230d4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -105,6 +105,7 @@ endmacro() if(OPENMC_USE_NCRYSTAL) find_package(NCrystal REQUIRED) + message(STATUS "Found NCrystal: ${NCrystal_DIR} (version ${NCrystal_VERSION})") endif() #=============================================================================== From a5b811ec899470f9b0bb3e387f0f672ae8490321 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Thu, 20 Oct 2022 13:25:05 +0200 Subject: [PATCH 10/42] Use new NCrystal getFlattenedComposition method to better support more complicated materials --- openmc/material.py | 63 +++++++++++++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 17 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 15d1f8ac9..de4a22cb8 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -335,15 +335,23 @@ 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. - + def from_ncrystal(cls, cfg, material_id=None, name=''): + """Create material from NCrystal configuration string. Density, + temperature, and material composition, and (ultimately) thermal neutron + scattering, will be automatically be provided by NCrystal based on this + string. The name and material_id parameters are simply passed on to the + Material constructor. + Parameters ---------- cfg : str NCrystal configuration string + material_id : int, optional + Unique identifier for the material. If not specified, an identifier will + automatically be assigned. + name : str, optional + Name of the material. If not specified, the name will be the empty + string. Returns ------- @@ -353,22 +361,43 @@ class Material(IDManagerMixin): """ import NCrystal - nc_mat = NCrystal.createInfo(cfg) - nc_comp = nc_mat.getComposition() + + def openmc_natabund( Z ): + #nc_mat.getFlattenedComposition might need natural abundancies. + #This call-back function is used so NCrystal can flatten composition + #using OpenMC's natural abundancies. In practice this function will + #only get invoked in the unlikely case where a material is specified + #by referring both to natural elements and specific isotopes of the + #same element. + elem_name = openmc.data.ATOMIC_SYMBOL.get( Z, None ) + if not elem_name: + raise ValueError( f'Element with Z={Z} is not known' ) + l = [] + for iso_name,abund in openmc.data.isotopes( elem_name ): + l.append( ( int(iso_name[ len(elem_name) : ]), abund ) ) + return l + + flat_compos = nc_mat.getFlattenedComposition( preferNaturalElements = True, + naturalAbundProvider = openmc_natabund ) # Create the Material - material = cls() + material = cls( material_id = material_id, + name = name, + temperature = nc_mat.getTemperature() ) - for frac, atom in nc_comp: - if not atom.isNaturalElement(): - raise ValueError('NCrystal-OpenMC interface only works with natural elements for now.') - material.add_element(atom.elementName(), frac, 'ao') + for Z, A_vals in flat_compos: + elemname = openmc.data.ATOMIC_SYMBOL.get(Z,None) + if not elemname: + raise ValueError(f'Element with Z={Z} is not known') + for A, frac in A_vals: + if A: + material.add_nuclide( elemname + str(A), frac, 'ao' ) + else: + material.add_element( elemname, frac, 'ao' ) - material._ncrystal_cfg = cfg - material._density_units = "g/cm3" - material._density = nc_mat.getDensity() - material.temperature = nc_mat.getTemperature() + material.set_density( 'g/cm3', nc_mat.getDensity() ) + material._ncrystal_cfg = NCrystal.normaliseCfg( cfg ) return material @@ -1186,7 +1215,7 @@ class Material(IDManagerMixin): 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.") + raise ValueError("NCrystal materials are not compatible with macroscopic cross sections.") element.set("cfg", str(self._ncrystal_cfg)) From f2e807e942ed421cd03fc87823c579e47dfcce3c Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Thu, 20 Oct 2022 13:26:45 +0200 Subject: [PATCH 11/42] Move NCrystalRNGWrapper data member into explicit private section as requested --- include/openmc/physics.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 23995e7af..806297c4e 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -111,8 +111,6 @@ void split_particle(Particle& p); //============================================================================== class NCrystalRNGWrapper : public NCrystal::RNGStream { - uint64_t* openmc_seed_; - public: constexpr NCrystalRNGWrapper(uint64_t* seed) noexcept : openmc_seed_(seed) {} @@ -122,6 +120,8 @@ protected: return std::max( std::numeric_limits::min(), prn(openmc_seed_)); } +private: + uint64_t* openmc_seed_; }; #endif From aca4fcb42bbfb45bd3ded0f9f1bbf8b7107aa632 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 21 Oct 2022 10:30:05 +0200 Subject: [PATCH 12/42] NCrystal cfgstr printout now mentions material ID --- src/material.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/material.cpp b/src/material.cpp index 9531e9f43..32f5f2e86 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -63,7 +63,7 @@ Material::Material(pugi::xml_node node) #ifdef NCRYSTAL if (check_for_node(node, "cfg")) { ncrystal_cfg_ = get_node_value(node, "cfg"); - write_message(5, "NCrystal config string: '{}'", ncrystal_cfg_); + write_message(5, "NCrystal config string for material #{}: '{}'", this->id(), ncrystal_cfg_); ncrystal_mat_ = NCrystal::FactImpl::createScatter(ncrystal_cfg_); } #endif From a3dd55d8bc4b28323b73e9a36c52f950db6ea1fb Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 21 Oct 2022 15:04:57 +0200 Subject: [PATCH 13/42] first stab at gha-install-ncrystal.sh --- tools/ci/gha-install-ncrystal.sh | 62 ++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100755 tools/ci/gha-install-ncrystal.sh diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh new file mode 100755 index 000000000..c18087f23 --- /dev/null +++ b/tools/ci/gha-install-ncrystal.sh @@ -0,0 +1,62 @@ +#!/bin/bash +set -ex +cd $HOME + +#Use the NCrystal develop branch (in the near future we can move this to master): +git clone https://github.com/mctools/ncrystal --branch develop --single-branch --depth 1 ncrystal_src + +SRC_DIR="$PWD/ncrystal_src" +BLD_DIR="$PWD/ncrystal_bld" +INST_DIR="$PWD/ncrystal_inst" +PYTHON=$(which python3) + +CPU_COUNT=1 + +mkdir "$BLD_DIR" +cd ncrystal_bld + +#cmake -Dstatic=on .. && make 2>/dev/null && sudo make install + +cmake \ + "${SRC_DIR}" \ + -DBUILD_SHARED_LIBS=ON \ + -DNCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE=ON \ + -DMODIFY_RPATH=OFF \ + -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_EXAMPLES=OFF \ + -DINSTALL_SETUPSH=OFF \ + -DEMBED_DATA=ON \ + -DINSTALL_DATA=OFF \ + -DNO_DIRECT_PYMODINST=ON \ + -DPython3_EXECUTABLE="$PYTHON" + +make -j${CPU_COUNT:-1} +make install + +#Note: There is no "make test" or "make ctest" functionality for NCrystal +# yet. If it appears in the future, we should add it here. + + +#The next stuff is not pretty, but it is needed as a temporary +#workaround until NCrystal add better support for system-wide +#installations in CMake: + +cat < ./findncrystallib.py +import pathlib +libname = pathlib.Path('./cfg_ncrystal_libname.txt').read_text().strip() +libs = list(pathlib.Path("/usr").glob("**/lib*/**/%s"%libname)) +assert len(libs)==1 +pathlib.Path('ncrystal_liblocation.txt').write_text(str(libs[0])) +EOF + +python3 ./findncrystallib.py + +TMPNCRYSTALLIBLOCATION=$(cat ./ncrystal_liblocation.txt) + +cat < ./ncrystal_pypkg/NCrystal/_nclibpath.py +#File autogenerated for working with systemwide install of NCrystal: +import pathlib +liblocation = pathlib.Path('${TMPNCRYSTALLIBLOCATION}'.strip()) +EOF + +$PYTHON -m pip install ./ncrystal_pypkg/ -vv From e64bf6a7ffb1b33f732e804bf10ff376d2cba700 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Thu, 10 Nov 2022 12:31:45 -0300 Subject: [PATCH 14/42] Add installation instructions for NCrystal --- docs/source/usersguide/install.rst | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 74be9ba0b..9e5ba1ffd 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -271,6 +271,16 @@ Prerequisites cmake -DOPENMC_USE_DAGMC=on -DCMAKE_PREFIX_PATH=/path/to/dagmc/installation .. + * NCrystal_ library for defining materials with enhanced thermal neutron transport + + Adding this option allows the creation of materials from NCrystal, which + replaces the scattering kernel treatment of ACE files with a modular, on-the-fly approach. + To use it `install `_ and + `initialize `_ + NCrystal and turn on the option in the CMake configuration step: + + cmake -DOPENMC_USE_NCRYSTAL=on .. + * libMesh_ mesh library framework for numerical simulations of partial differential equations This optional dependency enables support for unstructured mesh tally @@ -294,6 +304,7 @@ Prerequisites .. _MOAB: https://bitbucket.org/fathomteam/moab .. _libMesh: https://libmesh.github.io/ .. _libpng: http://www.libpng.org/pub/png/libpng.html +.. _NCrystal: https://github.com/mctools/ncrystal Obtaining the Source -------------------- @@ -362,6 +373,12 @@ OPENMC_USE_DAGMC should also be defined as `DAGMC_ROOT` in the CMake configuration command. (Default: off) +OPENMC_USE_NCRYSTAL + Turns on support for NCrystal materials. NCrystal must be + `installed `_ and + `initialized `_. + (Default: off) + OPENMC_USE_LIBMESH Enables the use of unstructured mesh tallies with libMesh_. (Default: off) From 1f3c5ed6bc4b29259d300a0d9b0fac0decf0eff6 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Thu, 10 Nov 2022 13:04:47 -0300 Subject: [PATCH 15/42] Add use instructions and examples for NCrystal --- docs/source/usersguide/materials.rst | 37 ++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 8b43f1a2a..ab5a8c5af 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -101,6 +101,42 @@ you would need to add hydrogen and oxygen to a material and then assign the .. _usersguide_naming: +------------------ +Adding NCrystal materials +------------------ + +Additional support for thermal scattering can be added by using NCrystal_. +The :meth:`Material.from_ncrystal` class method generates a material object from +an `NCrystal configuration string `_. +Temperature, material composition and density are passed from the configuration string +and the `NCMAT file `_ +that define the material, e.g: + +:: + + mat = openmc.Material.from_ncrystal('Al_sg225.ncmat') + +defines a material containing polycrystalline alumnium, + +:: + + mat = openmc.Material.from_ncrystal(""""Ge_sg227.ncmat;dcutoff=0.5;mos=40arcsec; + dir1=@crys_hkl:5,1,1@lab:0,0,1; + dir2=@crys_hkl:0,-1,1@lab:0,1,0""") + +defines an oriented germanium single crystal with 40 arcsec mosaicity. + +NCrystal only handles low energy neutron interactions. Other interaction are +provided by standard ACE files. NCrystal_ comes with a `predefined library +https://github.com/mctools/ncrystal/wiki/Data-library`_ but more materials can +be added by creating NCMAT files or on-the-fly in the configuration string. + +.. warning:: Currently, NCrystal_ materials cannot be modified after created. + Density, temperature and composition should be defined in the + configuration string or the NCMAT file. + +.. _NCrystal: https://github.com/mctools/ncrystal + ------------------ Naming Conventions ------------------ @@ -222,3 +258,4 @@ been generated, you can tell OpenMC to use this file either by setting materials.cross_sections = '/path/to/cross_sections.xml' .. _MCNP: https://mcnp.lanl.gov/ + From cf96ddf8bbedac4d72e1730e313d5b7a9cd94284 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Thu, 17 Nov 2022 13:47:03 +0100 Subject: [PATCH 16/42] Basic introduction to NCrystal Included references to the main papers and documentation --- docs/source/methods/cross_sections.rst | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 3396d7f25..c700b3f2d 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -178,6 +178,24 @@ been selected. There are three methods available: section data is loaded for a single temperature and is used in the unresolved resonance and fast energy ranges. +------------------ +NCrystal materials +------------------ + +As an alternative of the standard thermal scattering treatment using :math:`S(\alpha,\beta)` +tables, OpenMC allows to create materials using NCrystal_. In addition to the regular thermal elastic, +and thermal inelastic processes, NCrystal allows the generation of models for materials that +cannot currently included in ACE files such as oriented single crystals (see the `NCrystal paper`_), +and further extend the physics `using plugins`_. Thermal scattering kernels are generated on +the fly from dynamic and structural data, or loaded from :math:`S(\alpha,\beta)` tables converted +from ENDF6 evaluations. These kernels are sampled in a direct way using a fast `rejection algorithm`_ +that does not require previous processing. A `large library` of materials is already included in +the NCrystal distribution, and new materials can be easily defined from scratch in the `NCMAT format` +or `combining existing files`. + +The compositions of the materials defined in NCrystal are passed on to OpenMC all other reactions +except for thermal neutron scattering are handled by continuous energy ACE libraries. + ---------------- Multi-Group Data ---------------- @@ -279,3 +297,10 @@ or even isotropic scattering. .. _ENDF/B data: https://www.nndc.bnl.gov/endf-b8.0/ .. _Leppanen: https://doi.org/10.1016/j.anucene.2009.03.019 .. _algorithms: http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package +.. _NCrystal: https://github.com/mctools/ncrystal +.. _NCrystal paper: https://doi.org/10.1016/j.cpc.2019.07.015 +.. _using plugins: https://doi.org/10.1016/j.cpc.2021.108082 +.. _rejection algorithm: https://doi.org/10.1016/j.jcp.2018.11.043 +.. _large library: https://github.com/mctools/ncrystal/wiki/Data-library +.. _NCMAT format: https://github.com/mctools/ncrystal/wiki/NCMAT-format +.. _combining existing files: https://github.com/mctools/ncrystal/wiki/Announcement-Release3.0.0#2-multiphase-materials From 24572a71f22b588eb628e855618a4a4a27a63e4c Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Thu, 17 Nov 2022 16:33:48 +0100 Subject: [PATCH 17/42] Added simple NCrystal test --- tests/regression_tests/ncrystal/__init__.py | 0 .../regression_tests/ncrystal/inputs_true.dat | 45 +++++ .../ncrystal/results_true.dat | 181 ++++++++++++++++++ tests/regression_tests/ncrystal/test.py | 70 +++++++ 4 files changed, 296 insertions(+) create mode 100644 tests/regression_tests/ncrystal/__init__.py create mode 100644 tests/regression_tests/ncrystal/inputs_true.dat create mode 100644 tests/regression_tests/ncrystal/results_true.dat create mode 100644 tests/regression_tests/ncrystal/test.py diff --git a/tests/regression_tests/ncrystal/__init__.py b/tests/regression_tests/ncrystal/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/ncrystal/inputs_true.dat b/tests/regression_tests/ncrystal/inputs_true.dat new file mode 100644 index 000000000..7a3f20eb3 --- /dev/null +++ b/tests/regression_tests/ncrystal/inputs_true.dat @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + + + fixed source + 100000 + 10 + + + 0 0 -20 + + + + 0.012 1.0 + + + + + + + 1 + + + 0.0 0.017453292519943295 0.03490658503988659 0.05235987755982989 0.06981317007977318 0.08726646259971647 0.10471975511965978 0.12217304763960307 0.13962634015954636 0.15707963267948966 0.17453292519943295 0.19198621771937624 0.20943951023931956 0.22689280275926285 0.24434609527920614 0.2617993877991494 0.2792526803190927 0.29670597283903605 0.3141592653589793 0.33161255787892263 0.3490658503988659 0.3665191429188092 0.3839724354387525 0.4014257279586958 0.4188790204786391 0.4363323129985824 0.4537856055185257 0.47123889803846897 0.4886921905584123 0.5061454830783556 0.5235987755982988 0.5410520681182421 0.5585053606381855 0.5759586531581288 0.5934119456780721 0.6108652381980153 0.6283185307179586 0.6457718232379019 0.6632251157578453 0.6806784082777885 0.6981317007977318 0.7155849933176751 0.7330382858376184 0.7504915783575618 0.767944870877505 0.7853981633974483 0.8028514559173916 0.8203047484373349 0.8377580409572782 0.8552113334772214 0.8726646259971648 0.8901179185171081 0.9075712110370514 0.9250245035569946 0.9424777960769379 0.9599310885968813 0.9773843811168246 0.9948376736367679 1.0122909661567112 1.0297442586766545 1.0471975511965976 1.064650843716541 1.0821041362364843 1.0995574287564276 1.117010721276371 1.1344640137963142 1.1519173063162575 1.1693705988362009 1.1868238913561442 1.2042771838760873 1.2217304763960306 1.239183768915974 1.2566370614359172 1.2740903539558606 1.2915436464758039 1.3089969389957472 1.3264502315156905 1.3439035240356338 1.361356816555577 1.3788101090755203 1.3962634015954636 1.413716694115407 1.4311699866353502 1.4486232791552935 1.4660765716752369 1.4835298641951802 1.5009831567151235 1.5184364492350666 1.53588974175501 1.5533430342749532 1.5707963267948966 1.5882496193148399 1.6057029118347832 1.6231562043547265 1.6406094968746698 1.6580627893946132 1.6755160819145565 1.6929693744344996 1.710422666954443 1.7278759594743862 1.7453292519943295 1.7627825445142729 1.7802358370342162 1.7976891295541595 1.8151424220741028 1.8325957145940461 1.8500490071139892 1.8675022996339325 1.8849555921538759 1.9024088846738192 1.9198621771937625 1.9373154697137058 1.9547687622336491 1.9722220547535925 1.9896753472735358 2.007128639793479 2.0245819323134224 2.0420352248333655 2.059488517353309 2.076941809873252 2.0943951023931953 2.111848394913139 2.129301687433082 2.1467549799530254 2.1642082724729685 2.181661564992912 2.199114857512855 2.2165681500327987 2.234021442552742 2.251474735072685 2.2689280275926285 2.2863813201125716 2.303834612632515 2.321287905152458 2.3387411976724017 2.356194490192345 2.3736477827122884 2.3911010752322315 2.4085543677521746 2.426007660272118 2.443460952792061 2.4609142453120048 2.478367537831948 2.4958208303518914 2.5132741228718345 2.530727415391778 2.548180707911721 2.5656340004316642 2.5830872929516078 2.600540585471551 2.6179938779914944 2.6354471705114375 2.652900463031381 2.670353755551324 2.6878070480712677 2.705260340591211 2.722713633111154 2.7401669256310974 2.7576202181510405 2.775073510670984 2.792526803190927 2.8099800957108707 2.827433388230814 2.8448866807507573 2.8623399732707004 2.8797932657906435 2.897246558310587 2.91469985083053 2.9321531433504737 2.949606435870417 2.9670597283903604 2.9845130209103035 3.001966313430247 3.01941960595019 3.036872898470133 3.0543261909900767 3.07177948351002 3.0892327760299634 3.1066860685499065 3.12413936106985 3.141592653589793 + + + 1 + + + 1 2 3 + current + + diff --git a/tests/regression_tests/ncrystal/results_true.dat b/tests/regression_tests/ncrystal/results_true.dat new file mode 100644 index 000000000..3ef0306fd --- /dev/null +++ b/tests/regression_tests/ncrystal/results_true.dat @@ -0,0 +1,181 @@ + surface polar low [rad] polar high [rad] cellfrom nuclide score mean std. dev. +0 1 0.00e+00 1.75e-02 1 total current 9.82e-01 1.43e-04 +1 1 1.75e-02 3.49e-02 1 total current 0.00e+00 0.00e+00 +2 1 3.49e-02 5.24e-02 1 total current 1.00e-06 1.00e-06 +3 1 5.24e-02 6.98e-02 1 total current 0.00e+00 0.00e+00 +4 1 6.98e-02 8.73e-02 1 total current 0.00e+00 0.00e+00 +5 1 8.73e-02 1.05e-01 1 total current 1.00e-06 1.00e-06 +6 1 1.05e-01 1.22e-01 1 total current 0.00e+00 0.00e+00 +7 1 1.22e-01 1.40e-01 1 total current 1.00e-06 1.00e-06 +8 1 1.40e-01 1.57e-01 1 total current 1.00e-06 1.00e-06 +9 1 1.57e-01 1.75e-01 1 total current 1.00e-06 1.00e-06 +10 1 1.75e-01 1.92e-01 1 total current 0.00e+00 0.00e+00 +11 1 1.92e-01 2.09e-01 1 total current 2.00e-06 1.33e-06 +12 1 2.09e-01 2.27e-01 1 total current 0.00e+00 0.00e+00 +13 1 2.27e-01 2.44e-01 1 total current 1.00e-06 1.00e-06 +14 1 2.44e-01 2.62e-01 1 total current 1.00e-06 1.00e-06 +15 1 2.62e-01 2.79e-01 1 total current 2.00e-06 1.33e-06 +16 1 2.79e-01 2.97e-01 1 total current 0.00e+00 0.00e+00 +17 1 2.97e-01 3.14e-01 1 total current 2.00e-06 2.00e-06 +18 1 3.14e-01 3.32e-01 1 total current 1.00e-06 1.00e-06 +19 1 3.32e-01 3.49e-01 1 total current 1.00e-06 1.00e-06 +20 1 3.49e-01 3.67e-01 1 total current 2.00e-06 1.33e-06 +21 1 3.67e-01 3.84e-01 1 total current 0.00e+00 0.00e+00 +22 1 3.84e-01 4.01e-01 1 total current 2.00e-06 1.33e-06 +23 1 4.01e-01 4.19e-01 1 total current 2.00e-06 1.33e-06 +24 1 4.19e-01 4.36e-01 1 total current 1.00e-06 1.00e-06 +25 1 4.36e-01 4.54e-01 1 total current 2.00e-06 2.00e-06 +26 1 4.54e-01 4.71e-01 1 total current 5.00e-06 2.24e-06 +27 1 4.71e-01 4.89e-01 1 total current 4.00e-06 1.63e-06 +28 1 4.89e-01 5.06e-01 1 total current 3.00e-06 1.53e-06 +29 1 5.06e-01 5.24e-01 1 total current 3.00e-06 1.53e-06 +30 1 5.24e-01 5.41e-01 1 total current 3.00e-06 1.53e-06 +31 1 5.41e-01 5.59e-01 1 total current 7.00e-06 2.13e-06 +32 1 5.59e-01 5.76e-01 1 total current 3.00e-06 1.53e-06 +33 1 5.76e-01 5.93e-01 1 total current 3.00e-06 1.53e-06 +34 1 5.93e-01 6.11e-01 1 total current 2.00e-06 1.33e-06 +35 1 6.11e-01 6.28e-01 1 total current 2.00e-06 1.33e-06 +36 1 6.28e-01 6.46e-01 1 total current 3.00e-06 2.13e-06 +37 1 6.46e-01 6.63e-01 1 total current 3.00e-06 1.53e-06 +38 1 6.63e-01 6.81e-01 1 total current 2.00e-06 1.33e-06 +39 1 6.81e-01 6.98e-01 1 total current 3.00e-06 1.53e-06 +40 1 6.98e-01 7.16e-01 1 total current 1.00e-06 1.00e-06 +41 1 7.16e-01 7.33e-01 1 total current 6.00e-06 2.67e-06 +42 1 7.33e-01 7.50e-01 1 total current 6.00e-06 2.21e-06 +43 1 7.50e-01 7.68e-01 1 total current 7.00e-06 3.35e-06 +44 1 7.68e-01 7.85e-01 1 total current 6.00e-06 2.67e-06 +45 1 7.85e-01 8.03e-01 1 total current 6.00e-06 2.21e-06 +46 1 8.03e-01 8.20e-01 1 total current 7.00e-06 2.13e-06 +47 1 8.20e-01 8.38e-01 1 total current 5.00e-06 2.24e-06 +48 1 8.38e-01 8.55e-01 1 total current 3.00e-06 1.53e-06 +49 1 8.55e-01 8.73e-01 1 total current 8.00e-06 3.27e-06 +50 1 8.73e-01 8.90e-01 1 total current 7.00e-06 2.13e-06 +51 1 8.90e-01 9.08e-01 1 total current 6.00e-06 2.21e-06 +52 1 9.08e-01 9.25e-01 1 total current 4.00e-06 2.21e-06 +53 1 9.25e-01 9.42e-01 1 total current 1.20e-05 3.27e-06 +54 1 9.42e-01 9.60e-01 1 total current 8.00e-06 3.89e-06 +55 1 9.60e-01 9.77e-01 1 total current 1.20e-05 4.42e-06 +56 1 9.77e-01 9.95e-01 1 total current 7.00e-06 3.67e-06 +57 1 9.95e-01 1.01e+00 1 total current 1.20e-05 3.89e-06 +58 1 1.01e+00 1.03e+00 1 total current 9.00e-06 2.33e-06 +59 1 1.03e+00 1.05e+00 1 total current 6.00e-06 2.67e-06 +60 1 1.05e+00 1.06e+00 1 total current 7.00e-06 2.13e-06 +61 1 1.06e+00 1.08e+00 1 total current 5.00e-06 2.24e-06 +62 1 1.08e+00 1.10e+00 1 total current 1.20e-05 2.91e-06 +63 1 1.10e+00 1.12e+00 1 total current 1.30e-05 3.35e-06 +64 1 1.12e+00 1.13e+00 1 total current 9.00e-06 3.14e-06 +65 1 1.13e+00 1.15e+00 1 total current 1.20e-05 3.89e-06 +66 1 1.15e+00 1.17e+00 1 total current 6.00e-06 2.21e-06 +67 1 1.17e+00 1.19e+00 1 total current 5.11e-03 4.20e-05 +68 1 1.19e+00 1.20e+00 1 total current 8.00e-06 3.27e-06 +69 1 1.20e+00 1.22e+00 1 total current 1.30e-05 4.48e-06 +70 1 1.22e+00 1.24e+00 1 total current 1.20e-05 3.89e-06 +71 1 1.24e+00 1.26e+00 1 total current 1.50e-05 4.28e-06 +72 1 1.26e+00 1.27e+00 1 total current 7.00e-06 3.00e-06 +73 1 1.27e+00 1.29e+00 1 total current 1.60e-05 3.71e-06 +74 1 1.29e+00 1.31e+00 1 total current 9.00e-06 3.79e-06 +75 1 1.31e+00 1.33e+00 1 total current 1.30e-05 2.60e-06 +76 1 1.33e+00 1.34e+00 1 total current 1.60e-05 3.06e-06 +77 1 1.34e+00 1.36e+00 1 total current 1.40e-05 2.67e-06 +78 1 1.36e+00 1.38e+00 1 total current 1.40e-05 6.86e-06 +79 1 1.38e+00 1.40e+00 1 total current 1.50e-05 4.28e-06 +80 1 1.40e+00 1.41e+00 1 total current 3.26e-03 6.93e-05 +81 1 1.41e+00 1.43e+00 1 total current 1.40e-05 3.71e-06 +82 1 1.43e+00 1.45e+00 1 total current 1.30e-05 3.00e-06 +83 1 1.45e+00 1.47e+00 1 total current 1.30e-05 3.67e-06 +84 1 1.47e+00 1.48e+00 1 total current 1.70e-05 5.59e-06 +85 1 1.48e+00 1.50e+00 1 total current 1.70e-05 3.67e-06 +86 1 1.50e+00 1.52e+00 1 total current 1.40e-05 3.40e-06 +87 1 1.52e+00 1.54e+00 1 total current 2.50e-05 4.53e-06 +88 1 1.54e+00 1.55e+00 1 total current 1.30e-05 5.39e-06 +89 1 1.55e+00 1.57e+00 1 total current 1.80e-05 4.67e-06 +90 1 1.57e+00 1.59e+00 1 total current 1.40e-05 4.52e-06 +91 1 1.59e+00 1.61e+00 1 total current 1.70e-05 4.23e-06 +92 1 1.61e+00 1.62e+00 1 total current 1.40e-05 3.71e-06 +93 1 1.62e+00 1.64e+00 1 total current 1.00e-05 2.11e-06 +94 1 1.64e+00 1.66e+00 1 total current 2.00e-05 4.22e-06 +95 1 1.66e+00 1.68e+00 1 total current 2.30e-05 5.59e-06 +96 1 1.68e+00 1.69e+00 1 total current 1.70e-05 6.51e-06 +97 1 1.69e+00 1.71e+00 1 total current 1.30e-05 3.00e-06 +98 1 1.71e+00 1.73e+00 1 total current 1.50e-05 4.01e-06 +99 1 1.73e+00 1.75e+00 1 total current 1.70e-05 3.96e-06 +100 1 1.75e+00 1.76e+00 1 total current 1.80e-05 5.12e-06 +101 1 1.76e+00 1.78e+00 1 total current 2.50e-05 6.54e-06 +102 1 1.78e+00 1.80e+00 1 total current 1.80e-05 3.59e-06 +103 1 1.80e+00 1.82e+00 1 total current 1.50e-05 4.01e-06 +104 1 1.82e+00 1.83e+00 1 total current 1.10e-05 4.07e-06 +105 1 1.83e+00 1.85e+00 1 total current 1.50e-05 4.01e-06 +106 1 1.85e+00 1.87e+00 1 total current 1.90e-05 4.82e-06 +107 1 1.87e+00 1.88e+00 1 total current 2.20e-05 3.89e-06 +108 1 1.88e+00 1.90e+00 1 total current 2.10e-05 3.79e-06 +109 1 1.90e+00 1.92e+00 1 total current 1.50e-05 3.42e-06 +110 1 1.92e+00 1.94e+00 1 total current 2.20e-05 5.12e-06 +111 1 1.94e+00 1.95e+00 1 total current 2.10e-05 5.86e-06 +112 1 1.95e+00 1.97e+00 1 total current 2.60e-05 4.27e-06 +113 1 1.97e+00 1.99e+00 1 total current 2.20e-05 4.16e-06 +114 1 1.99e+00 2.01e+00 1 total current 2.40e-05 5.42e-06 +115 1 2.01e+00 2.02e+00 1 total current 1.60e-05 4.52e-06 +116 1 2.02e+00 2.04e+00 1 total current 1.30e-05 3.35e-06 +117 1 2.04e+00 2.06e+00 1 total current 1.90e-05 4.07e-06 +118 1 2.06e+00 2.08e+00 1 total current 1.30e-05 3.00e-06 +119 1 2.08e+00 2.09e+00 1 total current 1.40e-05 4.00e-06 +120 1 2.09e+00 2.11e+00 1 total current 3.10e-05 5.47e-06 +121 1 2.11e+00 2.13e+00 1 total current 2.30e-05 5.39e-06 +122 1 2.13e+00 2.15e+00 1 total current 2.20e-05 4.67e-06 +123 1 2.15e+00 2.16e+00 1 total current 1.70e-05 4.48e-06 +124 1 2.16e+00 2.18e+00 1 total current 1.60e-05 4.00e-06 +125 1 2.18e+00 2.20e+00 1 total current 1.80e-05 4.16e-06 +126 1 2.20e+00 2.22e+00 1 total current 1.80e-05 4.42e-06 +127 1 2.22e+00 2.23e+00 1 total current 1.80e-05 5.54e-06 +128 1 2.23e+00 2.25e+00 1 total current 1.90e-05 4.33e-06 +129 1 2.25e+00 2.27e+00 1 total current 1.00e-05 3.33e-06 +130 1 2.27e+00 2.29e+00 1 total current 1.80e-05 4.42e-06 +131 1 2.29e+00 2.30e+00 1 total current 4.15e-03 5.40e-05 +132 1 2.30e+00 2.32e+00 1 total current 1.90e-05 2.33e-06 +133 1 2.32e+00 2.34e+00 1 total current 2.30e-05 3.96e-06 +134 1 2.34e+00 2.36e+00 1 total current 1.90e-05 3.48e-06 +135 1 2.36e+00 2.37e+00 1 total current 1.60e-05 3.71e-06 +136 1 2.37e+00 2.39e+00 1 total current 1.60e-05 4.27e-06 +137 1 2.39e+00 2.41e+00 1 total current 2.30e-05 4.48e-06 +138 1 2.41e+00 2.43e+00 1 total current 2.10e-05 3.48e-06 +139 1 2.43e+00 2.44e+00 1 total current 1.60e-05 2.21e-06 +140 1 2.44e+00 2.46e+00 1 total current 9.00e-06 2.77e-06 +141 1 2.46e+00 2.48e+00 1 total current 1.30e-05 3.00e-06 +142 1 2.48e+00 2.50e+00 1 total current 2.70e-05 3.67e-06 +143 1 2.50e+00 2.51e+00 1 total current 1.90e-05 4.07e-06 +144 1 2.51e+00 2.53e+00 1 total current 1.20e-05 4.16e-06 +145 1 2.53e+00 2.55e+00 1 total current 1.30e-05 2.13e-06 +146 1 2.55e+00 2.57e+00 1 total current 1.10e-05 2.33e-06 +147 1 2.57e+00 2.58e+00 1 total current 1.50e-05 3.07e-06 +148 1 2.58e+00 2.60e+00 1 total current 1.20e-05 2.49e-06 +149 1 2.60e+00 2.62e+00 1 total current 1.80e-05 5.54e-06 +150 1 2.62e+00 2.64e+00 1 total current 1.30e-05 3.67e-06 +151 1 2.64e+00 2.65e+00 1 total current 1.60e-05 3.40e-06 +152 1 2.65e+00 2.67e+00 1 total current 7.00e-06 3.35e-06 +153 1 2.67e+00 2.69e+00 1 total current 1.00e-05 2.98e-06 +154 1 2.69e+00 2.71e+00 1 total current 7.00e-06 3.35e-06 +155 1 2.71e+00 2.72e+00 1 total current 1.20e-05 2.91e-06 +156 1 2.72e+00 2.74e+00 1 total current 9.00e-06 2.33e-06 +157 1 2.74e+00 2.76e+00 1 total current 1.00e-05 3.33e-06 +158 1 2.76e+00 2.78e+00 1 total current 1.10e-05 3.14e-06 +159 1 2.78e+00 2.79e+00 1 total current 1.00e-05 3.33e-06 +160 1 2.79e+00 2.81e+00 1 total current 1.40e-05 4.76e-06 +161 1 2.81e+00 2.83e+00 1 total current 8.00e-06 2.91e-06 +162 1 2.83e+00 2.84e+00 1 total current 5.00e-06 2.69e-06 +163 1 2.84e+00 2.86e+00 1 total current 6.00e-06 2.21e-06 +164 1 2.86e+00 2.88e+00 1 total current 5.00e-06 1.67e-06 +165 1 2.88e+00 2.90e+00 1 total current 4.00e-06 2.21e-06 +166 1 2.90e+00 2.91e+00 1 total current 7.00e-06 2.13e-06 +167 1 2.91e+00 2.93e+00 1 total current 6.00e-06 2.67e-06 +168 1 2.93e+00 2.95e+00 1 total current 7.00e-06 2.13e-06 +169 1 2.95e+00 2.97e+00 1 total current 5.00e-06 1.67e-06 +170 1 2.97e+00 2.98e+00 1 total current 3.00e-06 1.53e-06 +171 1 2.98e+00 3.00e+00 1 total current 6.00e-06 2.21e-06 +172 1 3.00e+00 3.02e+00 1 total current 3.00e-06 1.53e-06 +173 1 3.02e+00 3.04e+00 1 total current 1.00e-05 2.98e-06 +174 1 3.04e+00 3.05e+00 1 total current 2.00e-06 1.33e-06 +175 1 3.05e+00 3.07e+00 1 total current 1.00e-06 1.00e-06 +176 1 3.07e+00 3.09e+00 1 total current 2.00e-06 1.33e-06 +177 1 3.09e+00 3.11e+00 1 total current 0.00e+00 0.00e+00 +178 1 3.11e+00 3.12e+00 1 total current 1.00e-06 1.00e-06 +179 1 3.12e+00 3.14e+00 1 total current 0.00e+00 0.00e+00 \ No newline at end of file diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py new file mode 100644 index 000000000..51a0f07b7 --- /dev/null +++ b/tests/regression_tests/ncrystal/test.py @@ -0,0 +1,70 @@ +import numpy as np +import openmc +from tests.testing_harness import PyAPITestHarness + +def compute_angular_distribution(cfg, E0, N): + """Return a openmc.model.Model() object for a monoenergetic pencil + beam hitting a 1 mm sphere filled with the material defined by + the cfg string, and compute the angular distribution""" + + # Material definition + + m1 = openmc.Material.from_ncrystal(cfg) + materials = openmc.Materials([m1]) + + # Geometry definition + + sample_sphere = openmc.Sphere(r=0.1) + outer_sphere = openmc.Sphere(r=100, boundary_type="vacuum") + cell1 = openmc.Cell(region= -sample_sphere, fill=m1) + cell2_region= +sample_sphere&-outer_sphere + cell2 = openmc.Cell(region= cell2_region, fill=None) + uni1 = openmc.Universe(cells=[cell1, cell2]) + geometry = openmc.Geometry(uni1) + + # Source definition + + source = openmc.Source() + source.space = openmc.stats.Point(xyz = (0,0,-20)) + source.angle = openmc.stats.Monodirectional(reference_uvw = (0,0,1)) + source.energy = openmc.stats.Discrete([E0], [1.0]) + + # Execution settings + + settings = openmc.Settings() + settings.source = source + settings.run_mode = "fixed source" + settings.batches = 10 + settings.particles = N + + # Tally definition + + tally1 = openmc.Tally(name="angular distribution") + tally1.scores = ["current"] + filter1 = openmc.filter.SurfaceFilter(sample_sphere) + filter2 = openmc.filter.PolarFilter(np.linspace(0,np.pi,180+1)) + filter3 = openmc.filter.CellFromFilter(cell1) + tally1.filters = [filter1, filter2, filter3] + tallies = openmc.Tallies([tally1]) + + return openmc.model.Model(geometry, materials, settings, tallies) + + +class NCrystalTest(PyAPITestHarness): + def _get_results(self): + """Digest info in the statepoint and return as a string.""" + + # Read the statepoint file. + sp = openmc.StatePoint(self._sp_name) + tal = sp.get_tally(name='angular distribution') + df = tal.get_pandas_dataframe() + return df.to_string() + +def test_ncrystal(): + NParticles = 100000 + T = 293.6 # K + E0 = 0.012 # eV + cfg = 'Al_sg225.ncmat' + test = compute_angular_distribution(cfg, E0, NParticles) + harness = NCrystalTest('statepoint.10.h5', model=test) + harness.main() From 54745f507641489bc82acbfb4303a728c7bf4fca Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Fri, 18 Nov 2022 11:11:37 -0300 Subject: [PATCH 18/42] Changed example NCrystal cfg strings --- docs/source/usersguide/materials.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index ab5a8c5af..1886a772b 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -114,13 +114,13 @@ that define the material, e.g: :: - mat = openmc.Material.from_ncrystal('Al_sg225.ncmat') + mat = openmc.Material.from_ncrystal('Al_sg225.ncmat;temp=300K') defines a material containing polycrystalline alumnium, :: - mat = openmc.Material.from_ncrystal(""""Ge_sg227.ncmat;dcutoff=0.5;mos=40arcsec; + mat = openmc.Material.from_ncrystal("""Ge_sg227.ncmat;dcutoff=0.5;mos=40arcsec; dir1=@crys_hkl:5,1,1@lab:0,0,1; dir2=@crys_hkl:0,-1,1@lab:0,1,0""") From 91f939a391032db72220e4be715b888d0f8db287 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Fri, 18 Nov 2022 18:26:31 +0100 Subject: [PATCH 19/42] Add flag for NCrystal --- include/openmc/settings.h | 2 ++ openmc/lib/__init__.py | 3 +++ src/settings.cpp | 7 +++++++ tests/regression_tests/ncrystal/test.py | 7 +++++++ 4 files changed, 19 insertions(+) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 9bb447005..b01c75d45 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -20,6 +20,8 @@ namespace openmc { // Global variable declarations //============================================================================== +extern "C" const bool NCRYSTAL_ENABLED; + namespace settings { // Boolean flags diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index c14b0d9c2..55e2ec11b 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -42,6 +42,9 @@ else: def _dagmc_enabled(): return c_bool.in_dll(_dll, "DAGMC_ENABLED").value +def _ncrystal_enabled(): + return c_bool.in_dll(_dll, "NCRYSTAL_ENABLED").value + def _coord_levels(): return c_int.in_dll(_dll, "n_coord_levels").value diff --git a/src/settings.cpp b/src/settings.cpp index 87c24eedd..be7228d68 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -36,6 +36,13 @@ 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 diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py index 51a0f07b7..0a47d634d 100644 --- a/tests/regression_tests/ncrystal/test.py +++ b/tests/regression_tests/ncrystal/test.py @@ -1,7 +1,14 @@ import numpy as np import openmc +import openmc.lib +import pytest + from tests.testing_harness import PyAPITestHarness +pytestmark = pytest.mark.skipif( + not openmc.lib._ncrystal_enabled(), + reason="NCrystal materials are not enabled.") + def compute_angular_distribution(cfg, E0, N): """Return a openmc.model.Model() object for a monoenergetic pencil beam hitting a 1 mm sphere filled with the material defined by From 5f8e64b463aa2e0509d2b9955dc7de6686d5497b Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Fri, 18 Nov 2022 17:47:21 -0300 Subject: [PATCH 20/42] Add NCrystal vairables to ci.yml --- .github/workflows/ci.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e267c0d8e..321f18874 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,6 +29,7 @@ jobs: mpi: [n, y] omp: [n, y] dagmc: [n] + ncrystal: [n] libmesh: [n] event: [n] vectfit: [n] @@ -47,6 +48,10 @@ jobs: python-version: '3.10' mpi: y omp: y + - ncrystal: y + python-version: '3.10' + mpi: n + omp: n - libmesh: y python-version: '3.10' mpi: y @@ -64,7 +69,7 @@ jobs: omp: n mpi: y name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }}, - mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }}, + mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }}, ncrystal=${{ matrix.ncrystal }}, libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }} vectfit=${{ matrix.vectfit }})" @@ -73,6 +78,7 @@ jobs: PHDF5: ${{ matrix.mpi }} OMP: ${{ matrix.omp }} DAGMC: ${{ matrix.dagmc }} + NCRYSTAL: ${{ matrix.ncrystal }} EVENT: ${{ matrix.event }} VECTFIT: ${{ matrix.vectfit }} LIBMESH: ${{ matrix.libmesh }} From 87ed26babee2c7b69c229529e146d335a90fd952 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Fri, 18 Nov 2022 17:48:19 -0300 Subject: [PATCH 21/42] Add option to install NCrystal --- tools/ci/gha-install.sh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index aa40eb90b..e5390be1e 100755 --- a/tools/ci/gha-install.sh +++ b/tools/ci/gha-install.sh @@ -17,6 +17,11 @@ if [[ $DAGMC = 'y' ]]; then ./tools/ci/gha-install-dagmc.sh fi +# Install NCrystal if needed +if [[ $NCRYSTAL = 'y' ]]; then + ./tools/ci/gha-install-ncrystal.sh +fi + # Install vectfit for WMP generation if needed if [[ $VECTFIT = 'y' ]]; then ./tools/ci/gha-install-vectfit.sh From 8cf9dbe320baf4ff98eb57979f6e249e69ee0424 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Fri, 18 Nov 2022 22:28:04 +0100 Subject: [PATCH 22/42] Fixed in NCrystal docs --- docs/source/methods/cross_sections.rst | 4 ++-- docs/source/usersguide/materials.rst | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index c700b3f2d..64e4fac3d 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -190,8 +190,8 @@ and further extend the physics `using plugins`_. Thermal scattering kernels are the fly from dynamic and structural data, or loaded from :math:`S(\alpha,\beta)` tables converted from ENDF6 evaluations. These kernels are sampled in a direct way using a fast `rejection algorithm`_ that does not require previous processing. A `large library` of materials is already included in -the NCrystal distribution, and new materials can be easily defined from scratch in the `NCMAT format` -or `combining existing files`. +the NCrystal distribution, and new materials can be easily defined from scratch in the `NCMAT format`_ +or `combining existing files`_. The compositions of the materials defined in NCrystal are passed on to OpenMC all other reactions except for thermal neutron scattering are handled by continuous energy ACE libraries. diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 1886a772b..1eaae846c 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -127,8 +127,8 @@ defines a material containing polycrystalline alumnium, defines an oriented germanium single crystal with 40 arcsec mosaicity. NCrystal only handles low energy neutron interactions. Other interaction are -provided by standard ACE files. NCrystal_ comes with a `predefined library -https://github.com/mctools/ncrystal/wiki/Data-library`_ but more materials can +provided by standard ACE files. NCrystal_ comes with a `predefined library +`_ but more materials can be added by creating NCMAT files or on-the-fly in the configuration string. .. warning:: Currently, NCrystal_ materials cannot be modified after created. From b9e8f1324e72dc4730185b5f1f350eea62ab6757 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Fri, 18 Nov 2022 22:40:46 +0100 Subject: [PATCH 23/42] Attempt to fix NCrystal installation script --- tools/ci/gha-install-ncrystal.sh | 4 ++++ tools/ci/gha-install.py | 8 ++++++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index c18087f23..5f577b8da 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -28,6 +28,7 @@ cmake \ -DEMBED_DATA=ON \ -DINSTALL_DATA=OFF \ -DNO_DIRECT_PYMODINST=ON \ + -DCMAKE_INSTALL_PREFIX="${INST_DIR}" \ -DPython3_EXECUTABLE="$PYTHON" make -j${CPU_COUNT:-1} @@ -60,3 +61,6 @@ liblocation = pathlib.Path('${TMPNCRYSTALLIBLOCATION}'.strip()) EOF $PYTHON -m pip install ./ncrystal_pypkg/ -vv + +setup ${INST_DIR}/setup.sh + diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index 4c69ba70b..f3a29273f 100644 --- a/tools/ci/gha-install.py +++ b/tools/ci/gha-install.py @@ -19,7 +19,7 @@ def which(program): return None -def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False): +def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False, ncrystal=False): # Create build directory and change to it shutil.rmtree('build', ignore_errors=True) os.mkdir('build') @@ -54,6 +54,9 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False): libmesh_path = os.environ.get('HOME') + '/LIBMESH' cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + libmesh_path) + if ncrystal: + cmake_cmd.append('-DOPENMC_USE_NCRYSTAL=ON') + # Build in coverage mode for coverage testing cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on') @@ -70,10 +73,11 @@ def main(): mpi = (os.environ.get('MPI') == 'y') phdf5 = (os.environ.get('PHDF5') == 'y') dagmc = (os.environ.get('DAGMC') == 'y') + ncrystal = (os.environ.get('NCRYSTAL') == 'y') libmesh = (os.environ.get('LIBMESH') == 'y') # Build and install - install(omp, mpi, phdf5, dagmc, libmesh) + install(omp, mpi, phdf5, dagmc, libmesh, ncrystal) if __name__ == '__main__': main() From d7ab491f473e87f5865e365bd5f9fbcfcd24943e Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Sat, 19 Nov 2022 08:57:00 +0100 Subject: [PATCH 24/42] Try fix of NCrystal install --- tools/ci/gha-install-ncrystal.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 5f577b8da..78a1a44f2 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -45,7 +45,7 @@ make install cat < ./findncrystallib.py import pathlib libname = pathlib.Path('./cfg_ncrystal_libname.txt').read_text().strip() -libs = list(pathlib.Path("/usr").glob("**/lib*/**/%s"%libname)) +libs = list(pathlib.Path("${INST_DIR}").glob("**/lib*/**/%s"%libname)) assert len(libs)==1 pathlib.Path('ncrystal_liblocation.txt').write_text(str(libs[0])) EOF From b199db9cb9c6de05117c2749c21ce5ea3d474f7d Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Sat, 19 Nov 2022 09:06:42 +0100 Subject: [PATCH 25/42] Update gha-install-ncrystal.sh --- tools/ci/gha-install-ncrystal.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 78a1a44f2..715e4da94 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -62,5 +62,5 @@ EOF $PYTHON -m pip install ./ncrystal_pypkg/ -vv -setup ${INST_DIR}/setup.sh +source ${INST_DIR}/setup.sh From 39d0c5cd61224b2a80b8c6ee75597588e132842b Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Sat, 19 Nov 2022 09:38:34 +0100 Subject: [PATCH 26/42] Try to initialize with setup.sh --- tools/ci/gha-install-ncrystal.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 715e4da94..8e75cc60a 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -24,7 +24,7 @@ cmake \ -DMODIFY_RPATH=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DBUILD_EXAMPLES=OFF \ - -DINSTALL_SETUPSH=OFF \ + -DINSTALL_SETUPSH=ON \ -DEMBED_DATA=ON \ -DINSTALL_DATA=OFF \ -DNO_DIRECT_PYMODINST=ON \ From 9e1abb328429072d9349fd5c01ae6249c438d187 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Sat, 19 Nov 2022 10:47:26 +0100 Subject: [PATCH 27/42] Update gha-install-ncrystal.sh --- tools/ci/gha-install-ncrystal.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 8e75cc60a..c080c16d8 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -24,7 +24,7 @@ cmake \ -DMODIFY_RPATH=OFF \ -DCMAKE_BUILD_TYPE=Release \ -DBUILD_EXAMPLES=OFF \ - -DINSTALL_SETUPSH=ON \ + -DINSTALL_SETUPSH=OFF \ -DEMBED_DATA=ON \ -DINSTALL_DATA=OFF \ -DNO_DIRECT_PYMODINST=ON \ @@ -62,5 +62,5 @@ EOF $PYTHON -m pip install ./ncrystal_pypkg/ -vv -source ${INST_DIR}/setup.sh +# source ${INST_DIR}/setup.sh From 3e6bb9d5f669d2729fa96acdb4894ed1d338f2c5 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 08:51:23 -0300 Subject: [PATCH 28/42] Upgrade installation script to NCrystal 3.5.0 --- tools/ci/gha-install-ncrystal.sh | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index c080c16d8..18472dfa6 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -21,13 +21,12 @@ cmake \ "${SRC_DIR}" \ -DBUILD_SHARED_LIBS=ON \ -DNCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE=ON \ - -DMODIFY_RPATH=OFF \ + -DNCRYSTAL_MODIFY_RPATH=OFF \ -DCMAKE_BUILD_TYPE=Release \ - -DBUILD_EXAMPLES=OFF \ - -DINSTALL_SETUPSH=OFF \ - -DEMBED_DATA=ON \ - -DINSTALL_DATA=OFF \ - -DNO_DIRECT_PYMODINST=ON \ + -DNCRYSTAL_ENABLE_EXAMPLES=OFF \ + -DNCRYSTAL_ENABLE_SETUPSH=OFF \ + -DNCRYSTAL_ENABLE_DATA=OFF \ + -DNCRYSTAL_SKIP_PYMODINST=ON \ -DCMAKE_INSTALL_PREFIX="${INST_DIR}" \ -DPython3_EXECUTABLE="$PYTHON" @@ -62,5 +61,5 @@ EOF $PYTHON -m pip install ./ncrystal_pypkg/ -vv -# source ${INST_DIR}/setup.sh +eval $( "${INST_DIR}/bin/ncrystal-config --setup" ) From be330a11b6798850ddaa54635167c89030bfb511 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 09:15:45 -0300 Subject: [PATCH 29/42] Update gha-install-ncrystal.sh --- tools/ci/gha-install-ncrystal.sh | 31 ++++--------------------------- 1 file changed, 4 insertions(+), 27 deletions(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 18472dfa6..1bcc89627 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -25,8 +25,7 @@ cmake \ -DCMAKE_BUILD_TYPE=Release \ -DNCRYSTAL_ENABLE_EXAMPLES=OFF \ -DNCRYSTAL_ENABLE_SETUPSH=OFF \ - -DNCRYSTAL_ENABLE_DATA=OFF \ - -DNCRYSTAL_SKIP_PYMODINST=ON \ + -DNCRYSTAL_ENABLE_DATA=EMBED \ -DCMAKE_INSTALL_PREFIX="${INST_DIR}" \ -DPython3_EXECUTABLE="$PYTHON" @@ -36,30 +35,8 @@ make install #Note: There is no "make test" or "make ctest" functionality for NCrystal # yet. If it appears in the future, we should add it here. - -#The next stuff is not pretty, but it is needed as a temporary -#workaround until NCrystal add better support for system-wide -#installations in CMake: - -cat < ./findncrystallib.py -import pathlib -libname = pathlib.Path('./cfg_ncrystal_libname.txt').read_text().strip() -libs = list(pathlib.Path("${INST_DIR}").glob("**/lib*/**/%s"%libname)) -assert len(libs)==1 -pathlib.Path('ncrystal_liblocation.txt').write_text(str(libs[0])) -EOF - -python3 ./findncrystallib.py - -TMPNCRYSTALLIBLOCATION=$(cat ./ncrystal_liblocation.txt) - -cat < ./ncrystal_pypkg/NCrystal/_nclibpath.py -#File autogenerated for working with systemwide install of NCrystal: -import pathlib -liblocation = pathlib.Path('${TMPNCRYSTALLIBLOCATION}'.strip()) -EOF - -$PYTHON -m pip install ./ncrystal_pypkg/ -vv - eval $( "${INST_DIR}/bin/ncrystal-config --setup" ) +# Check installation worked + +ncrystal-config --setup From fff545ac79d39b6aad0ee4dd055ac1da38cc344a Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 09:27:12 -0300 Subject: [PATCH 30/42] Update gha-install-ncrystal.sh --- tools/ci/gha-install-ncrystal.sh | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 1bcc89627..390b7dba5 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -35,8 +35,14 @@ make install #Note: There is no "make test" or "make ctest" functionality for NCrystal # yet. If it appears in the future, we should add it here. -eval $( "${INST_DIR}/bin/ncrystal-config --setup" ) +# Output the configuration to the log + +"${INST_DIR}/bin/ncrystal-config" --setup + +# Change environmental variables + +eval $( "${INST_DIR}/bin/ncrystal-config" --setup ) # Check installation worked -ncrystal-config --setup +nctool --test From 8598ff7e4311c80c3dc7993a34eb6fc40976d7ed Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 09:45:27 -0300 Subject: [PATCH 31/42] Pass CMAKE_PREFIX_PATH from NCrystall installation --- tools/ci/gha-install.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index f3a29273f..95fe70abe 100644 --- a/tools/ci/gha-install.py +++ b/tools/ci/gha-install.py @@ -56,6 +56,8 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False, ncrys if ncrystal: cmake_cmd.append('-DOPENMC_USE_NCRYSTAL=ON') + ncrystal_path = os.environ.get('CMAKE_PREFIX_PATH') + cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + ncrystal_path) # Build in coverage mode for coverage testing cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on') From c98a12171010e00fb78572c72f7327b0de740718 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 09:56:10 -0300 Subject: [PATCH 32/42] Update gha-install.py --- tools/ci/gha-install.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index 95fe70abe..f5779c935 100644 --- a/tools/ci/gha-install.py +++ b/tools/ci/gha-install.py @@ -56,8 +56,8 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False, ncrys if ncrystal: cmake_cmd.append('-DOPENMC_USE_NCRYSTAL=ON') - ncrystal_path = os.environ.get('CMAKE_PREFIX_PATH') - cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + ncrystal_path) + ncrystal_cmake_path = os.environ.get('HOME')+'/ncrystal_inst/lib/cmake' + cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + ncrystal_cmake_path) # Build in coverage mode for coverage testing cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on') From 19bf624b34806efcffc52a8f18ad915624b8811b Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 11:03:41 -0300 Subject: [PATCH 33/42] Add NCrystal installation test --- tools/ci/gha-script.sh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tools/ci/gha-script.sh b/tools/ci/gha-script.sh index c791167e9..4c125b529 100755 --- a/tools/ci/gha-script.sh +++ b/tools/ci/gha-script.sh @@ -14,5 +14,10 @@ if [[ $EVENT == 'y' ]]; then args="${args} --event " fi +# Check NCrystal installation +if [[ $NCRYSTAL = 'y' ]]; then + nctool --test +fi + # Run regression and unit tests pytest --cov=openmc -v $args tests From 4391d601fde4f74f91d594238e98564632c695b8 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 11:06:06 -0300 Subject: [PATCH 34/42] find NCrystal --- cmake/OpenMCConfig.cmake.in | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index d0e2beb82..162e3aad7 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -9,6 +9,11 @@ if(@OPENMC_USE_DAGMC@) find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@) endif() +if(@OPENMC_USE_NCRYSTAL@) + find_package(NCrystal REQUIRED) + message(STATUS "Found NCrystal: ${NCrystal_DIR} (version ${NCrystal_VERSION})") +endif() + if(@OPENMC_USE_LIBMESH@) include(FindPkgConfig) list(APPEND CMAKE_PREFIX_PATH @LIBMESH_PREFIX@) From e626447a4e8dbce5ff8211b0ada661c329e3cd3d Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com> Date: Mon, 28 Nov 2022 11:18:32 -0300 Subject: [PATCH 35/42] Update gha-script.sh --- tools/ci/gha-script.sh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/ci/gha-script.sh b/tools/ci/gha-script.sh index 4c125b529..85fdef49f 100755 --- a/tools/ci/gha-script.sh +++ b/tools/ci/gha-script.sh @@ -16,6 +16,8 @@ fi # Check NCrystal installation if [[ $NCRYSTAL = 'y' ]]; then + # Change environmental variables + eval $( "${HOME}/ncrystal_inst/bin/ncrystal-config" --setup ) nctool --test fi From e7b4285bd5214d0028700c07897fe60919ef8689 Mon Sep 17 00:00:00 2001 From: Jose Ignacio Marquez Damian Date: Tue, 29 Nov 2022 11:21:04 +0100 Subject: [PATCH 36/42] Cleanup --- docs/source/methods/cross_sections.rst | 2 +- docs/source/usersguide/materials.rst | 2 +- tools/ci/gha-install-ncrystal.sh | 2 -- 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 64e4fac3d..0b87913e5 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -189,7 +189,7 @@ cannot currently included in ACE files such as oriented single crystals (see the and further extend the physics `using plugins`_. Thermal scattering kernels are generated on the fly from dynamic and structural data, or loaded from :math:`S(\alpha,\beta)` tables converted from ENDF6 evaluations. These kernels are sampled in a direct way using a fast `rejection algorithm`_ -that does not require previous processing. A `large library` of materials is already included in +that does not require previous processing. A `large library`_ of materials is already included in the NCrystal distribution, and new materials can be easily defined from scratch in the `NCMAT format`_ or `combining existing files`_. diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 1eaae846c..a3c05f421 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -106,7 +106,7 @@ Adding NCrystal materials ------------------ Additional support for thermal scattering can be added by using NCrystal_. -The :meth:`Material.from_ncrystal` class method generates a material object from +The :meth:`Material.from_ncrystal` class method generates a :class:`openmc.Material` object from an `NCrystal configuration string `_. Temperature, material composition and density are passed from the configuration string and the `NCMAT file `_ diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh index 390b7dba5..16f77e13e 100755 --- a/tools/ci/gha-install-ncrystal.sh +++ b/tools/ci/gha-install-ncrystal.sh @@ -15,8 +15,6 @@ CPU_COUNT=1 mkdir "$BLD_DIR" cd ncrystal_bld -#cmake -Dstatic=on .. && make 2>/dev/null && sudo make install - cmake \ "${SRC_DIR}" \ -DBUILD_SHARED_LIBS=ON \ From f2c28706a6eaf2206da5ff73867c279946b55f4b Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Wed, 7 Dec 2022 10:12:23 +0100 Subject: [PATCH 37/42] Apply suggestions from code review Co-authored-by: Paul Romano --- CMakeLists.txt | 2 -- docs/source/usersguide/materials.rst | 14 ++++---- openmc/material.py | 43 +++++++++++++------------ tests/regression_tests/ncrystal/test.py | 29 ++++++++--------- tools/ci/gha-install.py | 4 +-- tools/ci/gha-script.sh | 6 ++-- 6 files changed, 47 insertions(+), 51 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2815230d4..5c4d489e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -497,8 +497,6 @@ if(OPENMC_USE_NCRYSTAL) target_link_libraries(libopenmc NCrystal::NCrystal) endif() - - #=============================================================================== # openmc executable #=============================================================================== diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index a3c05f421..fbf9effd8 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -101,19 +101,17 @@ you would need to add hydrogen and oxygen to a material and then assign the .. _usersguide_naming: ------------------- +------------------------- Adding NCrystal materials ------------------- +------------------------- Additional support for thermal scattering can be added by using NCrystal_. The :meth:`Material.from_ncrystal` class method generates a :class:`openmc.Material` object from an `NCrystal configuration string `_. -Temperature, material composition and density are passed from the configuration string +Temperature, material composition, and density are passed from the configuration string and the `NCMAT file `_ -that define the material, e.g: +that define the material, e.g.:: -:: - mat = openmc.Material.from_ncrystal('Al_sg225.ncmat;temp=300K') defines a material containing polycrystalline alumnium, @@ -126,12 +124,12 @@ defines a material containing polycrystalline alumnium, defines an oriented germanium single crystal with 40 arcsec mosaicity. -NCrystal only handles low energy neutron interactions. Other interaction are +NCrystal only handles low energy neutron interactions. Other interactions are provided by standard ACE files. NCrystal_ comes with a `predefined library `_ but more materials can be added by creating NCMAT files or on-the-fly in the configuration string. -.. warning:: Currently, NCrystal_ materials cannot be modified after created. +.. warning:: Currently, NCrystal_ materials cannot be modified after they are created. Density, temperature and composition should be defined in the configuration string or the NCMAT file. diff --git a/openmc/material.py b/openmc/material.py index de4a22cb8..298731c92 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -336,11 +336,12 @@ class Material(IDManagerMixin): @classmethod def from_ncrystal(cls, cfg, material_id=None, name=''): - """Create material from NCrystal configuration string. Density, - temperature, and material composition, and (ultimately) thermal neutron - scattering, will be automatically be provided by NCrystal based on this - string. The name and material_id parameters are simply passed on to the - Material constructor. + """Create material from NCrystal configuration string. + + Density, temperature, and material composition, and (ultimately) thermal + neutron scattering will be automatically be provided by NCrystal based + on this string. The name and material_id parameters are simply passed on + to the Material constructor. Parameters ---------- @@ -363,41 +364,41 @@ class Material(IDManagerMixin): import NCrystal nc_mat = NCrystal.createInfo(cfg) - def openmc_natabund( Z ): + def openmc_natabund(Z): #nc_mat.getFlattenedComposition might need natural abundancies. #This call-back function is used so NCrystal can flatten composition #using OpenMC's natural abundancies. In practice this function will #only get invoked in the unlikely case where a material is specified #by referring both to natural elements and specific isotopes of the #same element. - elem_name = openmc.data.ATOMIC_SYMBOL.get( Z, None ) + elem_name = openmc.data.ATOMIC_SYMBOL.get(Z) if not elem_name: raise ValueError( f'Element with Z={Z} is not known' ) - l = [] - for iso_name,abund in openmc.data.isotopes( elem_name ): - l.append( ( int(iso_name[ len(elem_name) : ]), abund ) ) - return l + return [ + (int(iso_name[len(elem_name):]), abund) + for iso_name, abund in openmc.data.isotopes(elem_name) + ] - flat_compos = nc_mat.getFlattenedComposition( preferNaturalElements = True, - naturalAbundProvider = openmc_natabund ) + flat_compos = nc_mat.getFlattenedComposition(preferNaturalElements = True, + naturalAbundProvider=openmc_natabund) # Create the Material - material = cls( material_id = material_id, - name = name, - temperature = nc_mat.getTemperature() ) + material = cls(material_id=material_id, + name=name, + temperature=nc_mat.getTemperature()) for Z, A_vals in flat_compos: - elemname = openmc.data.ATOMIC_SYMBOL.get(Z,None) + elemname = openmc.data.ATOMIC_SYMBOL.get(Z) if not elemname: raise ValueError(f'Element with Z={Z} is not known') for A, frac in A_vals: if A: - material.add_nuclide( elemname + str(A), frac, 'ao' ) + material.add_nuclide(f'{elemname}{A}', frac) else: - material.add_element( elemname, frac, 'ao' ) + material.add_element(elemname, frac) - material.set_density( 'g/cm3', nc_mat.getDensity() ) - material._ncrystal_cfg = NCrystal.normaliseCfg( cfg ) + material.set_density('g/cm3', nc_mat.getDensity()) + material._ncrystal_cfg = NCrystal.normaliseCfg(cfg) return material diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py index 0a47d634d..3137515d4 100644 --- a/tests/regression_tests/ncrystal/test.py +++ b/tests/regression_tests/ncrystal/test.py @@ -10,7 +10,7 @@ pytestmark = pytest.mark.skipif( reason="NCrystal materials are not enabled.") def compute_angular_distribution(cfg, E0, N): - """Return a openmc.model.Model() object for a monoenergetic pencil + """Return an openmc.Model() object for a monoenergetic pencil beam hitting a 1 mm sphere filled with the material defined by the cfg string, and compute the angular distribution""" @@ -23,17 +23,16 @@ def compute_angular_distribution(cfg, E0, N): sample_sphere = openmc.Sphere(r=0.1) outer_sphere = openmc.Sphere(r=100, boundary_type="vacuum") - cell1 = openmc.Cell(region= -sample_sphere, fill=m1) - cell2_region= +sample_sphere&-outer_sphere - cell2 = openmc.Cell(region= cell2_region, fill=None) - uni1 = openmc.Universe(cells=[cell1, cell2]) - geometry = openmc.Geometry(uni1) + cell1 = openmc.Cell(region=-sample_sphere, fill=m1) + cell2_region= +sample_sphere & -outer_sphere + cell2 = openmc.Cell(region=cell2_region, fill=None) + geometry = openmc.Geometry([cell1, cell2]) # Source definition source = openmc.Source() - source.space = openmc.stats.Point(xyz = (0,0,-20)) - source.angle = openmc.stats.Monodirectional(reference_uvw = (0,0,1)) + source.space = openmc.stats.Point((0, 0, -20)) + source.angle = openmc.stats.Monodirectional(reference_uvw=(0, 0, 1)) source.energy = openmc.stats.Discrete([E0], [1.0]) # Execution settings @@ -48,13 +47,13 @@ def compute_angular_distribution(cfg, E0, N): tally1 = openmc.Tally(name="angular distribution") tally1.scores = ["current"] - filter1 = openmc.filter.SurfaceFilter(sample_sphere) - filter2 = openmc.filter.PolarFilter(np.linspace(0,np.pi,180+1)) - filter3 = openmc.filter.CellFromFilter(cell1) + filter1 = openmc.SurfaceFilter(sample_sphere) + filter2 = openmc.PolarFilter(np.linspace(0, np.pi, 180+1)) + filter3 = openmc.CellFromFilter(cell1) tally1.filters = [filter1, filter2, filter3] tallies = openmc.Tallies([tally1]) - return openmc.model.Model(geometry, materials, settings, tallies) + return openmc.Model(geometry, materials, settings, tallies) class NCrystalTest(PyAPITestHarness): @@ -62,9 +61,9 @@ class NCrystalTest(PyAPITestHarness): """Digest info in the statepoint and return as a string.""" # Read the statepoint file. - sp = openmc.StatePoint(self._sp_name) - tal = sp.get_tally(name='angular distribution') - df = tal.get_pandas_dataframe() + with openmc.StatePoint(self._sp_name) as sp: + tal = sp.get_tally(name='angular distribution') + df = tal.get_pandas_dataframe() return df.to_string() def test_ncrystal(): diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index f5779c935..ac10bb3f6 100644 --- a/tools/ci/gha-install.py +++ b/tools/ci/gha-install.py @@ -56,8 +56,8 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False, ncrys if ncrystal: cmake_cmd.append('-DOPENMC_USE_NCRYSTAL=ON') - ncrystal_cmake_path = os.environ.get('HOME')+'/ncrystal_inst/lib/cmake' - cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + ncrystal_cmake_path) + ncrystal_cmake_path = os.environ.get('HOME') + '/ncrystal_inst/lib/cmake' + cmake_cmd.append(f'-DCMAKE_PREFIX_PATH={ncrystal_cmake_path}') # Build in coverage mode for coverage testing cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on') diff --git a/tools/ci/gha-script.sh b/tools/ci/gha-script.sh index 85fdef49f..1f1c3a1eb 100755 --- a/tools/ci/gha-script.sh +++ b/tools/ci/gha-script.sh @@ -16,9 +16,9 @@ fi # Check NCrystal installation if [[ $NCRYSTAL = 'y' ]]; then - # Change environmental variables - eval $( "${HOME}/ncrystal_inst/bin/ncrystal-config" --setup ) - nctool --test + # Change environmental variables + eval $( "${HOME}/ncrystal_inst/bin/ncrystal-config" --setup ) + nctool --test fi # Run regression and unit tests From a86e74f59a27236a4705ab0b4dc8c79736cd2c56 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Wed, 7 Dec 2022 12:51:31 +0100 Subject: [PATCH 38/42] Rename compute_angular_distribution to pencil_beam_model --- tests/regression_tests/ncrystal/test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py index 3137515d4..2a789dac7 100644 --- a/tests/regression_tests/ncrystal/test.py +++ b/tests/regression_tests/ncrystal/test.py @@ -9,7 +9,7 @@ pytestmark = pytest.mark.skipif( not openmc.lib._ncrystal_enabled(), reason="NCrystal materials are not enabled.") -def compute_angular_distribution(cfg, E0, N): +def pencil_beam_model(cfg, E0, N): """Return an openmc.Model() object for a monoenergetic pencil beam hitting a 1 mm sphere filled with the material defined by the cfg string, and compute the angular distribution""" @@ -71,6 +71,6 @@ def test_ncrystal(): T = 293.6 # K E0 = 0.012 # eV cfg = 'Al_sg225.ncmat' - test = compute_angular_distribution(cfg, E0, NParticles) + test = pencil_beam_model(cfg, E0, NParticles) harness = NCrystalTest('statepoint.10.h5', model=test) harness.main() From 62af77310897b6de1ef43906a781da0c9fcfe15f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Jan 2023 18:16:45 +0700 Subject: [PATCH 39/42] Refactor NCrystal interface code to ncrystal_interface.h/cpp --- CMakeLists.txt | 1 + include/openmc/material.h | 25 ++----- include/openmc/ncrystal_interface.h | 91 +++++++++++++++++++++++ include/openmc/physics.h | 28 +------- include/openmc/settings.h | 8 +-- src/material.cpp | 39 +++------- src/ncrystal_interface.cpp | 108 ++++++++++++++++++++++++++++ src/physics.cpp | 30 ++------ src/settings.cpp | 12 +--- 9 files changed, 229 insertions(+), 113 deletions(-) create mode 100644 include/openmc/ncrystal_interface.h create mode 100644 src/ncrystal_interface.cpp 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 //============================================================================== From 4eda693133a3a0f6c55eea084ca7ddeb3ab0769d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Jan 2023 10:39:22 +0700 Subject: [PATCH 40/42] Style and spacing changes --- docs/source/methods/cross_sections.rst | 27 ++++++++++++++----------- docs/source/usersguide/install.rst | 11 +++++----- docs/source/usersguide/materials.rst | 20 +++++++++--------- include/openmc/physics.h | 5 +++-- tests/regression_tests/ncrystal/test.py | 16 +++++++++------ 5 files changed, 45 insertions(+), 34 deletions(-) diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 0b87913e5..2cafa8691 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -182,19 +182,22 @@ been selected. There are three methods available: NCrystal materials ------------------ -As an alternative of the standard thermal scattering treatment using :math:`S(\alpha,\beta)` -tables, OpenMC allows to create materials using NCrystal_. In addition to the regular thermal elastic, -and thermal inelastic processes, NCrystal allows the generation of models for materials that -cannot currently included in ACE files such as oriented single crystals (see the `NCrystal paper`_), -and further extend the physics `using plugins`_. Thermal scattering kernels are generated on -the fly from dynamic and structural data, or loaded from :math:`S(\alpha,\beta)` tables converted -from ENDF6 evaluations. These kernels are sampled in a direct way using a fast `rejection algorithm`_ -that does not require previous processing. A `large library`_ of materials is already included in -the NCrystal distribution, and new materials can be easily defined from scratch in the `NCMAT format`_ -or `combining existing files`_. +As an alternative of the standard thermal scattering treatment using +:math:`S(\alpha,\beta)` tables, OpenMC allows to create materials using +NCrystal_. In addition to the regular thermal elastic, and thermal inelastic +processes, NCrystal allows the generation of models for materials that cannot +currently included in ACE files such as oriented single crystals (see the +`NCrystal paper`_), and further extend the physics `using plugins`_. Thermal +scattering kernels are generated on the fly from dynamic and structural data, or +loaded from :math:`S(\alpha,\beta)` tables converted from ENDF6 evaluations. +These kernels are sampled in a direct way using a fast `rejection algorithm`_ +that does not require previous processing. A `large library`_ of materials is +already included in the NCrystal distribution, and new materials can be easily +defined from scratch in the `NCMAT format`_ or `combining existing files`_. -The compositions of the materials defined in NCrystal are passed on to OpenMC all other reactions -except for thermal neutron scattering are handled by continuous energy ACE libraries. +The compositions of the materials defined in NCrystal are passed on to OpenMC +all other reactions except for thermal neutron scattering are handled by +continuous energy ACE libraries. ---------------- Multi-Group Data diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 9e5ba1ffd..537dda824 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -274,9 +274,10 @@ Prerequisites * NCrystal_ library for defining materials with enhanced thermal neutron transport Adding this option allows the creation of materials from NCrystal, which - replaces the scattering kernel treatment of ACE files with a modular, on-the-fly approach. - To use it `install `_ and - `initialize `_ + replaces the scattering kernel treatment of ACE files with a modular, + on-the-fly approach. To use it `install + `_ and `initialize + `_ NCrystal and turn on the option in the CMake configuration step: cmake -DOPENMC_USE_NCRYSTAL=on .. @@ -374,8 +375,8 @@ OPENMC_USE_DAGMC (Default: off) OPENMC_USE_NCRYSTAL - Turns on support for NCrystal materials. NCrystal must be - `installed `_ and + Turns on support for NCrystal materials. NCrystal must be + `installed `_ and `initialized `_. (Default: off) diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index fbf9effd8..83af55805 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -105,12 +105,14 @@ you would need to add hydrogen and oxygen to a material and then assign the Adding NCrystal materials ------------------------- -Additional support for thermal scattering can be added by using NCrystal_. -The :meth:`Material.from_ncrystal` class method generates a :class:`openmc.Material` object from -an `NCrystal configuration string `_. -Temperature, material composition, and density are passed from the configuration string -and the `NCMAT file `_ -that define the material, e.g.:: +Additional support for thermal scattering can be added by using NCrystal_. The +:meth:`Material.from_ncrystal` class method generates a :class:`openmc.Material` +object from an `NCrystal configuration string +`_. +Temperature, material composition, and density are passed from the configuration +string and the `NCMAT file +`_ that define the +material, e.g.:: mat = openmc.Material.from_ncrystal('Al_sg225.ncmat;temp=300K') @@ -125,12 +127,12 @@ defines a material containing polycrystalline alumnium, defines an oriented germanium single crystal with 40 arcsec mosaicity. NCrystal only handles low energy neutron interactions. Other interactions are -provided by standard ACE files. NCrystal_ comes with a `predefined library +provided by standard ACE files. NCrystal_ comes with a `predefined library `_ but more materials can be added by creating NCMAT files or on-the-fly in the configuration string. -.. warning:: Currently, NCrystal_ materials cannot be modified after they are created. - Density, temperature and composition should be defined in the +.. warning:: Currently, NCrystal_ materials cannot be modified after they are created. + Density, temperature and composition should be defined in the configuration string or the NCMAT file. .. _NCrystal: https://github.com/mctools/ncrystal diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 382f0f0dd..6e7327d38 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -96,8 +96,9 @@ 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); diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py index 2a789dac7..cb6421b03 100644 --- a/tests/regression_tests/ncrystal/test.py +++ b/tests/regression_tests/ncrystal/test.py @@ -1,3 +1,5 @@ +from math import pi + import numpy as np import openmc import openmc.lib @@ -9,6 +11,7 @@ pytestmark = pytest.mark.skipif( not openmc.lib._ncrystal_enabled(), reason="NCrystal materials are not enabled.") + def pencil_beam_model(cfg, E0, N): """Return an openmc.Model() object for a monoenergetic pencil beam hitting a 1 mm sphere filled with the material defined by @@ -24,7 +27,7 @@ def pencil_beam_model(cfg, E0, N): sample_sphere = openmc.Sphere(r=0.1) outer_sphere = openmc.Sphere(r=100, boundary_type="vacuum") cell1 = openmc.Cell(region=-sample_sphere, fill=m1) - cell2_region= +sample_sphere & -outer_sphere + cell2_region = +sample_sphere & -outer_sphere cell2 = openmc.Cell(region=cell2_region, fill=None) geometry = openmc.Geometry([cell1, cell2]) @@ -48,7 +51,7 @@ def pencil_beam_model(cfg, E0, N): tally1 = openmc.Tally(name="angular distribution") tally1.scores = ["current"] filter1 = openmc.SurfaceFilter(sample_sphere) - filter2 = openmc.PolarFilter(np.linspace(0, np.pi, 180+1)) + filter2 = openmc.PolarFilter(np.linspace(0, pi, 180+1)) filter3 = openmc.CellFromFilter(cell1) tally1.filters = [filter1, filter2, filter3] tallies = openmc.Tallies([tally1]) @@ -66,11 +69,12 @@ class NCrystalTest(PyAPITestHarness): df = tal.get_pandas_dataframe() return df.to_string() + def test_ncrystal(): - NParticles = 100000 - T = 293.6 # K - E0 = 0.012 # eV + n_particles = 100000 + T = 293.6 # K + E0 = 0.012 # eV cfg = 'Al_sg225.ncmat' - test = pencil_beam_model(cfg, E0, NParticles) + test = pencil_beam_model(cfg, E0, n_particles) harness = NCrystalTest('statepoint.10.h5', model=test) harness.main() From c4227f7dc1b30879dd7492730e8dbb02b14006d5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Jan 2023 10:39:45 +0700 Subject: [PATCH 41/42] Change Material.from_ncrystal to accept kwargs --- openmc/material.py | 51 +++++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 26 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 68ec199c3..7844ec131 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -99,6 +99,10 @@ class Material(IDManagerMixin): [decay/sec]. .. versionadded:: 0.13.2 + ncrystal_cfg : str + NCrystal configuration string + + .. versionadded:: 0.13.3 """ @@ -141,7 +145,8 @@ class Material(IDManagerMixin): string += '{: <16}\n'.format('\tS(a,b) Tables') - string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg) + if self._ncrystal_cfg: + string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg) for sab in self._sab: string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab) @@ -222,6 +227,10 @@ class Material(IDManagerMixin): def volume(self): return self._volume + @property + def ncrystal_cfg(self): + return self._ncrystal_cfg + @name.setter def name(self, name: Optional[str]): if name is not None: @@ -335,9 +344,9 @@ class Material(IDManagerMixin): return material @classmethod - def from_ncrystal(cls, cfg, material_id=None, name=''): + def from_ncrystal(cls, cfg, **kwargs): """Create material from NCrystal configuration string. - + Density, temperature, and material composition, and (ultimately) thermal neutron scattering will be automatically be provided by NCrystal based on this string. The name and material_id parameters are simply passed on @@ -347,12 +356,8 @@ class Material(IDManagerMixin): ---------- cfg : str NCrystal configuration string - material_id : int, optional - Unique identifier for the material. If not specified, an identifier will - automatically be assigned. - name : str, optional - Name of the material. If not specified, the name will be the empty - string. + **kwargs + Keyword arguments passed to :class:`openmc.Material` Returns ------- @@ -371,26 +376,20 @@ class Material(IDManagerMixin): #only get invoked in the unlikely case where a material is specified #by referring both to natural elements and specific isotopes of the #same element. - elem_name = openmc.data.ATOMIC_SYMBOL.get(Z) - if not elem_name: - raise ValueError( f'Element with Z={Z} is not known' ) + elem_name = openmc.data.ATOMIC_SYMBOL[Z] return [ (int(iso_name[len(elem_name):]), abund) for iso_name, abund in openmc.data.isotopes(elem_name) ] - flat_compos = nc_mat.getFlattenedComposition(preferNaturalElements = True, - naturalAbundProvider=openmc_natabund) + flat_compos = nc_mat.getFlattenedComposition( + preferNaturalElements=True, naturalAbundProvider=openmc_natabund) # Create the Material - material = cls(material_id=material_id, - name=name, - temperature=nc_mat.getTemperature()) + material = cls(temperature=nc_mat.getTemperature(), **kwargs) for Z, A_vals in flat_compos: - elemname = openmc.data.ATOMIC_SYMBOL.get(Z) - if not elemname: - raise ValueError(f'Element with Z={Z} is not known') + elemname = openmc.data.ATOMIC_SYMBOL[Z] for A, frac in A_vals: if A: material.add_nuclide(f'{elemname}{A}', frac) @@ -1067,7 +1066,7 @@ class Material(IDManagerMixin): activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier return activity if by_nuclide else sum(activity.values()) - + def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False): """Returns the decay heat of the material or for each nuclide in the material in units of [W], [W/g] or [W/cm3]. @@ -1101,16 +1100,16 @@ class Material(IDManagerMixin): multiplier = 1 elif units == 'W/g': multiplier = 1.0 / self.get_mass_density() - - decayheat = {} + + decayheat = {} for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items(): decay_erg = openmc.data.decay_energy(nuclide) inv_seconds = openmc.data.decay_constant(nuclide) decay_erg *= openmc.data.JOULE_PER_EV decayheat[nuclide] = inv_seconds * decay_erg * 1e24 * atoms_per_bcm * multiplier - return decayheat if by_nuclide else sum(decayheat.values()) - + return decayheat if by_nuclide else sum(decayheat.values()) + def get_nuclide_atoms(self): """Return number of atoms of each nuclide in the material @@ -1651,4 +1650,4 @@ class Materials(cv.CheckedList): tree = ET.parse(path) root = tree.getroot() - return cls.from_xml_element(root) \ No newline at end of file + return cls.from_xml_element(root) From e83ced88dc58d0a9fcf83f04390a5284b30d06b9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Jan 2023 01:17:53 +0700 Subject: [PATCH 42/42] Change ncrystal_max_energy to be a global constant --- include/openmc/ncrystal_interface.h | 3 +++ include/openmc/settings.h | 2 -- src/material.cpp | 2 +- src/physics.cpp | 2 +- src/settings.cpp | 1 - 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/openmc/ncrystal_interface.h b/include/openmc/ncrystal_interface.h index 36038d50c..5a3882df9 100644 --- a/include/openmc/ncrystal_interface.h +++ b/include/openmc/ncrystal_interface.h @@ -20,6 +20,9 @@ namespace openmc { extern "C" const bool NCRYSTAL_ENABLED; +//! Energy in [eV] to switch between NCrystal and ENDF +constexpr double NCRYSTAL_MAX_ENERGY {5.0}; + //============================================================================== // Wrapper class an NCrystal material //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 9cd430deb..90f1e8c29 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -91,8 +91,6 @@ 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 diff --git a/src/material.cpp b/src/material.cpp index 2d55d6a04..74a8d7e35 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -801,7 +801,7 @@ void Material::calculate_neutron_xs(Particle& p) const // Calculate NCrystal cross section double ncrystal_xs = -1.0; - if (ncrystal_mat_ && p.E() < settings::ncrystal_max_energy) { + if (ncrystal_mat_ && p.E() < NCRYSTAL_MAX_ENERGY) { ncrystal_xs = ncrystal_mat_.xs(p); } diff --git a/src/physics.cpp b/src/physics.cpp index 46bf37aa9..7cb8040f6 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -142,7 +142,7 @@ void sample_neutron_reaction(Particle& p) // 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) { + if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { ncrystal_mat.scatter(p); } else { scatter(p, i_nuclide); diff --git a/src/settings.cpp b/src/settings.cpp index 9174d4114..2fda352ce 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -99,7 +99,6 @@ 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};