Refactor NCrystal interface code to ncrystal_interface.h/cpp

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

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<const NCrystal::ProcImpl::Process> 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<int> nuclide_; //!< Indices in nuclides vector
-  vector<int> element_; //!< Indices in elements vector
-#ifdef NCRYSTAL
-  std::string ncrystal_cfg_; //!< NCrystal configuration string
-  std::shared_ptr<const NCrystal::ProcImpl::Process> ncrystal_mat_;
-#endif
+  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]

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 <cstdint> // for uint64_t
+#include <limits>  // for numeric_limits
+#include <string>
+
+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<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

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<double>(
-      std::numeric_limits<double>::min(), prn(openmc_seed_));
-  }
-private:
-  uint64_t* openmc_seed_;
-};
-#endif
-
 } // namespace openmc
 
 #endif // OPENMC_PHYSICS_H

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
 
 //==============================================================================

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
     }
 
     // ======================================================================

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<double>(
+      std::numeric_limits<double>::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

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;

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
 
 //==============================================================================