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 }} diff --git a/CMakeLists.txt b/CMakeLists.txt index 0e29c71cb..4737ff47b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF) option(OPENMC_USE_MPI "Enable MPI" OFF) option(OPENMC_USE_MCPL "Enable MCPL" 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") @@ -99,6 +100,15 @@ macro(find_package_write_status pkg) endif() endmacro() +#=============================================================================== +# NCrystal Scattering Support +#=============================================================================== + +if(OPENMC_USE_NCRYSTAL) + find_package(NCrystal REQUIRED) + message(STATUS "Found NCrystal: ${NCrystal_DIR} (version ${NCrystal_VERSION})") +endif() + #=============================================================================== # DAGMC Geometry Support - need DAGMC/MOAB #=============================================================================== @@ -338,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 @@ -498,6 +509,11 @@ if (OPENMC_USE_MCPL) target_link_libraries(libopenmc MCPL::mcpl) endif() +if(OPENMC_USE_NCRYSTAL) + target_compile_definitions(libopenmc PRIVATE NCRYSTAL) + target_link_libraries(libopenmc NCrystal::NCrystal) +endif() + #=============================================================================== # Log build info that this executable can report later #=============================================================================== diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index 1e3508305..756fe26dc 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@) diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 3396d7f25..2cafa8691 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -178,6 +178,27 @@ 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 +300,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 diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 74be9ba0b..537dda824 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -271,6 +271,17 @@ 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 +305,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 +374,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) diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 8b43f1a2a..83af55805 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 :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') + +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 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 they are 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/ + diff --git a/include/openmc/file_utils.h b/include/openmc/file_utils.h index cb52b33c8..9612896be 100644 --- a/include/openmc/file_utils.h +++ b/include/openmc/file_utils.h @@ -14,7 +14,8 @@ namespace openmc { inline bool dir_exists(const std::string& path) { struct stat s; - if (stat(path.c_str(), &s) != 0) return false; + if (stat(path.c_str(), &s) != 0) + return false; return s.st_mode & S_IFDIR; } diff --git a/include/openmc/material.h b/include/openmc/material.h index 81f4e1421..8db9bb628 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -12,6 +12,7 @@ #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" @@ -151,12 +152,17 @@ public: //! \return Temperature in [K] double temperature() const; + //! Get pointer to NCrystal material object + //! \return Pointer to NCrystal material object + 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 + 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..5a3882df9 --- /dev/null +++ b/include/openmc/ncrystal_interface.h @@ -0,0 +1,94 @@ +#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; + +//! Energy in [eV] to switch between NCrystal and ENDF +constexpr double NCRYSTAL_MAX_ENERGY {5.0}; + +//============================================================================== +// 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 262b3a884..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/openmc/lib/__init__.py b/openmc/lib/__init__.py index 1e337fc07..eef18be88 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/openmc/material.py b/openmc/material.py index aa778cab6..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 """ @@ -118,6 +122,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 +145,9 @@ class Material(IDManagerMixin): string += '{: <16}\n'.format('\tS(a,b) Tables') + 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) @@ -219,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: @@ -331,6 +343,64 @@ class Material(IDManagerMixin): return material + @classmethod + 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 + to the Material constructor. + + Parameters + ---------- + cfg : str + NCrystal configuration string + **kwargs + Keyword arguments passed to :class:`openmc.Material` + + Returns + ------- + openmc.Material + Material instance + + """ + + import NCrystal + nc_mat = NCrystal.createInfo(cfg) + + 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[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) + + # Create the Material + material = cls(temperature=nc_mat.getTemperature(), **kwargs) + + for Z, A_vals in flat_compos: + elemname = openmc.data.ATOMIC_SYMBOL[Z] + for A, frac in A_vals: + if A: + material.add_nuclide(f'{elemname}{A}', frac) + else: + material.add_element(elemname, frac) + + material.set_density('g/cm3', nc_mat.getDensity()) + material._ncrystal_cfg = NCrystal.normaliseCfg(cfg) + + return material + def add_volume_information(self, volume_calc): """Add volume information to a material. @@ -405,6 +475,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 +682,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() @@ -990,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]. @@ -1024,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 @@ -1182,6 +1258,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 with 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)) @@ -1566,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) diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 2094fd551..a7bd86095 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -91,7 +91,8 @@ Library::Library(pugi::xml_node node, const std::string& directory) // Non-member functions //============================================================================== -void read_cross_sections_xml() { +void read_cross_sections_xml() +{ pugi::xml_document doc; std::string filename = settings::path_input + "materials.xml"; // Check if materials.xml exists diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 83e8494de..eabdc9f70 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -40,7 +40,8 @@ void update_universe_cell_count(int32_t a, int32_t b) } } -void read_geometry_xml() { +void read_geometry_xml() +{ // Display output message write_message("Reading geometry XML file...", 5); diff --git a/src/initialize.cpp b/src/initialize.cpp index 0c137795b..1f9e40548 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -106,7 +106,8 @@ int openmc_init(int argc, char* argv[], const void* intracomm) openmc::openmc_set_seed(DEFAULT_SEED); // Read XML input files - if (!read_model_xml()) read_separate_xml_files(); + if (!read_model_xml()) + read_separate_xml_files(); // Write some initial output under the header if needed initial_output(); @@ -299,7 +300,8 @@ int parse_command_line(int argc, char* argv[]) return 0; } -bool read_model_xml() { +bool read_model_xml() +{ std::string model_filename = settings::path_input.empty() ? "." : settings::path_input; @@ -314,7 +316,8 @@ bool read_model_xml() { model_filename += "/model.xml"; // if this file doesn't exist, stop here - if (!file_exists(model_filename)) return false; + if (!file_exists(model_filename)) + return false; // try to process the path input as an XML file pugi::xml_document doc; @@ -327,7 +330,7 @@ bool read_model_xml() { // Read settings if (!check_for_node(root, "settings")) { - fatal_error("No node present in the model.xml file."); + fatal_error("No node present in the model.xml file."); } auto settings_root = root.child("settings"); @@ -343,13 +346,15 @@ bool read_model_xml() { title(); } - write_message(fmt::format("Reading model XML file '{}' ...", model_filename), 5); + write_message( + fmt::format("Reading model XML file '{}' ...", model_filename), 5); read_settings_xml(settings_root); // If other XML files are present, display warning // that they will be ignored - auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml", "tallies.xml", "plots.xml"}; + auto other_inputs = {"materials.xml", "geometry.xml", "settings.xml", + "tallies.xml", "plots.xml"}; for (const auto& input : other_inputs) { if (file_exists(settings::path_input + input)) { warning((fmt::format("Other XML file input(s) are present. These files " @@ -384,7 +389,7 @@ bool read_model_xml() { if (check_for_node(root, "tallies")) read_tallies_xml(root.child("tallies")); - // Initialize distribcell_filters + // Initialize distribcell_filters prepare_distribcell(); if (check_for_node(root, "plots")) @@ -416,7 +421,8 @@ void read_separate_xml_files() read_plots_xml(); } -void initial_output() { +void initial_output() +{ // write initial output if (settings::run_mode == RunMode::PLOTTING) { // Read plots.xml if it exists diff --git a/src/material.cpp b/src/material.cpp index f2e23aede..74a8d7e35 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -60,6 +60,13 @@ Material::Material(pugi::xml_node node) name_ = get_node_value(node, "name"); } + if (check_for_node(node, "cfg")) { + auto cfg = get_node_value(node, "cfg"); + write_message( + 5, "NCrystal config string for material #{}: '{}'", this->id(), cfg); + ncrystal_mat_ = NCrystalMat(cfg); + } + if (check_for_node(node, "depletable")) { depletable_ = get_node_value_bool(node, "depletable"); } @@ -367,8 +374,8 @@ void Material::finalize() this->init_thermal(); } -// Normalize density -this->normalize_density(); + // Normalize density + this->normalize_density(); } void Material::normalize_density() @@ -792,6 +799,12 @@ void Material::calculate_neutron_xs(Particle& p) const // Initialize position in i_sab_nuclides int j = 0; + // Calculate NCrystal cross section + double ncrystal_xs = -1.0; + if (ncrystal_mat_ && p.E() < NCRYSTAL_MAX_ENERGY) { + ncrystal_xs = ncrystal_mat_.xs(p); + } + // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { // ====================================================================== @@ -830,10 +843,16 @@ 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)}; + 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); + + // If NCrystal is being used, update micro cross section cache + if (ncrystal_xs >= 0.0) { + data::nuclides[i_nuclide]->calculate_elastic_xs(p); + ncrystal_update_micro(ncrystal_xs, micro); + } } // ====================================================================== @@ -1264,7 +1283,8 @@ double density_effect(const vector& f, const vector& e_b_sq, return delta - w_sq * (1.0 - beta_sq); } -void read_materials_xml() { +void read_materials_xml() +{ write_message("Reading materials XML file...", 5); pugi::xml_document doc; 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 7d391fdf6..7cb8040f6 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" @@ -140,7 +141,12 @@ void sample_neutron_reaction(Particle& p) // Sample a scattering reaction and determine the secondary energy of the // exiting neutron - scatter(p, i_nuclide); + const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); + if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { + ncrystal_mat.scatter(p); + } else { + scatter(p, i_nuclide); + } // Advance URR seed stream 'N' times after energy changes if (p.E() != p.E_last()) { diff --git a/src/plot.cpp b/src/plot.cpp index 1a62fe263..fac13a1f4 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -124,7 +124,8 @@ extern "C" int openmc_plot_geometry() return 0; } -void read_plots_xml() { +void read_plots_xml() +{ // Check if plots.xml exists; this is only necessary when the plot runmode is // initiated. Otherwise, we want to read plots.xml because it may be called // later via the API. In that case, its ok for a plots.xml to not exist diff --git a/src/settings.cpp b/src/settings.cpp index 6386ce52b..2fda352ce 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -214,7 +214,8 @@ void get_run_parameters(pugi::xml_node node_base) } } -void read_settings_xml() { +void read_settings_xml() +{ using namespace settings; using namespace pugi; // Check if settings.xml exists diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 1c1f274d7..94bfc6cd4 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -705,7 +705,8 @@ std::string Tally::nuclide_name(int nuclide_idx) const // Non-member functions //============================================================================== -void read_tallies_xml() { +void read_tallies_xml() +{ // Check if tallies.xml exists. If not, just return since it is optional std::string filename = settings::path_input + "tallies.xml"; if (!file_exists(filename)) diff --git a/src/thermal.cpp b/src/thermal.cpp index 982bf226a..741f89ed1 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -119,22 +119,24 @@ ThermalScattering::ThermalScattering( if (!found) { // If no pairs found, check if the desired temperature falls within // bounds' tolerance - if (std::abs(T - temps_available[0]) <= settings::temperature_tolerance){ - if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[0])) == - temps_to_read.end()) { - temps_to_read.push_back(std::round(temps_available[0])); - }} - else if (std::abs(T - temps_available[n - 1]) <= settings::temperature_tolerance){ - if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[n - 1])) == - temps_to_read.end()){ - temps_to_read.push_back(std::round(temps_available[n - 1])); - }} - else { + if (std::abs(T - temps_available[0]) <= + settings::temperature_tolerance) { + if (std::find(temps_to_read.begin(), temps_to_read.end(), + std::round(temps_available[0])) == temps_to_read.end()) { + temps_to_read.push_back(std::round(temps_available[0])); + } + } else if (std::abs(T - temps_available[n - 1]) <= + settings::temperature_tolerance) { + if (std::find(temps_to_read.begin(), temps_to_read.end(), + std::round(temps_available[n - 1])) == temps_to_read.end()) { + temps_to_read.push_back(std::round(temps_available[n - 1])); + } + } else { fatal_error( - fmt::format("Nuclear data library does not contain cross " - "sections for {} at temperatures that bound {} K.", - name_, std::round(T))); - } + fmt::format("Nuclear data library does not contain cross " + "sections for {} at temperatures that bound {} K.", + name_, std::round(T))); + } } } } 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..cb6421b03 --- /dev/null +++ b/tests/regression_tests/ncrystal/test.py @@ -0,0 +1,80 @@ +from math import pi + +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 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""" + + # 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) + geometry = openmc.Geometry([cell1, cell2]) + + # Source definition + + source = openmc.Source() + 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 + + 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.SurfaceFilter(sample_sphere) + filter2 = openmc.PolarFilter(np.linspace(0, pi, 180+1)) + filter3 = openmc.CellFromFilter(cell1) + tally1.filters = [filter1, filter2, filter3] + tallies = openmc.Tallies([tally1]) + + return openmc.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. + 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(): + n_particles = 100000 + T = 293.6 # K + E0 = 0.012 # eV + cfg = 'Al_sg225.ncmat' + test = pencil_beam_model(cfg, E0, n_particles) + harness = NCrystalTest('statepoint.10.h5', model=test) + harness.main() diff --git a/tools/ci/gha-install-ncrystal.sh b/tools/ci/gha-install-ncrystal.sh new file mode 100755 index 000000000..16f77e13e --- /dev/null +++ b/tools/ci/gha-install-ncrystal.sh @@ -0,0 +1,46 @@ +#!/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 \ + "${SRC_DIR}" \ + -DBUILD_SHARED_LIBS=ON \ + -DNCRYSTAL_NOTOUCH_CMAKE_BUILD_TYPE=ON \ + -DNCRYSTAL_MODIFY_RPATH=OFF \ + -DCMAKE_BUILD_TYPE=Release \ + -DNCRYSTAL_ENABLE_EXAMPLES=OFF \ + -DNCRYSTAL_ENABLE_SETUPSH=OFF \ + -DNCRYSTAL_ENABLE_DATA=EMBED \ + -DCMAKE_INSTALL_PREFIX="${INST_DIR}" \ + -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. + +# 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 + +nctool --test diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index 83e5afac7..f4b2fbb31 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,11 @@ 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') + 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') @@ -70,10 +75,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() diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index 75fa0ca70..4042977c9 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 diff --git a/tools/ci/gha-script.sh b/tools/ci/gha-script.sh index c791167e9..1f1c3a1eb 100755 --- a/tools/ci/gha-script.sh +++ b/tools/ci/gha-script.sh @@ -14,5 +14,12 @@ if [[ $EVENT == 'y' ]]; then args="${args} --event " fi +# Check NCrystal installation +if [[ $NCRYSTAL = 'y' ]]; then + # Change environmental variables + eval $( "${HOME}/ncrystal_inst/bin/ncrystal-config" --setup ) + nctool --test +fi + # Run regression and unit tests pytest --cov=openmc -v $args tests