Added hooks for NCrystal

paulromano · Jan 12, 2023 · Sep 13, 2021 · Sep 13, 2022 · Sep 13, 2022 · Sep 13, 2022
commit 0e0cd0a22d059159a8799906f57d12cec99c1a89
diff --git 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
@@ -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 @@ class Material {
   std::string name_;                    //!< Name of material
   vector<int> nuclide_;                 //!< Indices in nuclides vector
   vector<int> element_;                 //!< Indices in elements vector
+#ifdef NCRYSTAL
+  std::string cfg_; //!< NCrystal configuration string
+  std::shared_ptr<const NCrystal::ProcImpl::Process> m_NCrystal_mat_;
+#endif
   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/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
@@ -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<double>( std::numeric_limits<double>::min(), prn(openmc_seed_) );
+  }
+
+};
+#endif
+
 } // namespace openmc
 
 #endif // OPENMC_PHYSICS_H
diff --git 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
@@ -118,6 +118,7 @@ def __init__(self, material_id=None, name='', temperature=None):
         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 @@ def __repr__(self):
 
         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 @@ def from_hdf5(cls, group: h5py.Group):
 
         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 @@ def add_nuclide(self, nuclide: str, percent: float, percent_type: str = 'ao'):
                   '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 @@ def add_element(self, element: str, percent: float, percent_type: str = 'ao',
             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 @@ def to_xml_element(self):
         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
@@ -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
@@ -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
@@ -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
 
 //==============================================================================