paulromano · Jan 12, 2023 · Sep 13, 2021 · Sep 13, 2022 · Sep 13, 2022 · Sep 13, 2022
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
+      <https://github.com/mctools/ncrystal/wiki/Get-NCrystal>`_ and `initialize
+      <https://github.com/mctools/ncrystal/wiki/Using-NCrystal#setting-up>`_
+      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 <https://github.com/mctools/ncrystal/wiki/Get-NCrystal>`_ and
+  `initialized <https://github.com/mctools/ncrystal/wiki/Using-NCrystal#setting-up>`_.
+  (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
@@ -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
+<https://github.com/mctools/ncrystal/wiki/Using-NCrystal#uniform-material-configuration-syntax>`_.
+Temperature, material composition, and density are passed from the configuration
+string and the `NCMAT file
+<https://github.com/mctools/ncrystal/wiki/NCMAT-format>`_ 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
+<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 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
@@ -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
@@ -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 @@ class Material {
   //! \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<int> nuclide_;                 //!< Indices in nuclides vector
   vector<int> element_;                 //!< Indices in elements vector
+  NCrystalMat ncrystal_mat_;            //!< NCrystal material object
   xt::xtensor<double, 1> 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
@@ -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 <cstdint> // for uint64_t
+#include <limits>  // for numeric_limits
+#include <string>
+
+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<const NCrystal::ProcImpl::Process>
+    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
@@ -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);
 

@@ -42,6 +42,9 @@
 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