Converted DagMC threadprivate variables to particle owned as well. Not quite as bad as I was expecting.

include/openmc/cell.h

@@ -110,7 +110,7 @@ public:
 
   //! Find the oncoming boundary of this cell.
   virtual std::pair<double, int32_t>
-  distance(Position r, Direction u, int32_t on_surface) const = 0;
+  distance(Position r, Direction u, int32_t on_surface, Particle* p) const = 0;
 
   //! Write all information needed to reconstruct the cell to an HDF5 group.
   //! \param group_id An HDF5 group id.
@@ -201,7 +201,7 @@ public:
   contains(Position r, Direction u, int32_t on_surface) const;
 
   std::pair<double, int32_t>
-  distance(Position r, Direction u, int32_t on_surface) const;
+  distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
 
   void to_hdf5(hid_t group_id) const;
 
@@ -245,7 +245,7 @@ public:
   bool contains(Position r, Direction u, int32_t on_surface) const;
 
   std::pair<double, int32_t>
-  distance(Position r, Direction u, int32_t on_surface) const;
+  distance(Position r, Direction u, int32_t on_surface, Particle* p) const;
 
   BoundingBox bounding_box() const;
 

include/openmc/dagmc.h

@@ -14,14 +14,6 @@ extern "C" const bool dagmc_enabled;
 
 namespace openmc {
 
-namespace simulation {
-
-extern moab::DagMC::RayHistory history; //!< facet history for DagMC particles
-extern Direction last_dir; //!< last direction passed to DagMC's ray_fire
-#pragma omp threadprivate(history, last_dir)
-
-}
-
 namespace model {
   extern moab::DagMC* DAG;
 }

include/openmc/particle.h

@@ -16,6 +16,10 @@
 #include "openmc/random_lcg.h"
 #include "openmc/tallies/filter_match.h"
 
+#ifdef DAGMC
+#include "DagMC.hpp"
+#endif
+
 
 namespace openmc {
 
@@ -349,6 +353,12 @@ public:
   double collision_distance_; // distance to particle's next closest collision
 
   int n_event_ {0}; // number of events executed in this particle's history
+
+  // DagMC state variables
+  #ifdef DAGMC
+  moab::DagMC::RayHistory history_;
+  Direction last_dir_;
+  #endif
 };
 
 } // namespace openmc

include/openmc/surface.h

@@ -11,6 +11,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/constants.h"
+#include "openmc/particle.h"
 #include "openmc/position.h"
 #include "dagmc.h"
 
@@ -113,8 +114,9 @@ public:
   //! Determine the direction of a ray reflected from the surface.
   //! \param[in] r The point at which the ray is incident.
   //! \param[in] u Incident direction of the ray
+  //! \param[inout] p Pointer to the particle
   //! \return Outgoing direction of the ray
-  virtual Direction reflect(Position r, Direction u) const;
+  virtual Direction reflect(Position r, Direction u, Particle* p) const;
 
   virtual Direction diffuse_reflect(Position r, Direction u,
     uint64_t* seed) const;
@@ -171,7 +173,7 @@ public:
   double evaluate(Position r) const;
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
-  Direction reflect(Position r, Direction u) const;
+  Direction reflect(Position r, Direction u, Particle* p) const;
 
   void to_hdf5(hid_t group_id) const;
 

src/cell.cpp

@@ -495,7 +495,7 @@ CSGCell::contains(Position r, Direction u, int32_t on_surface) const
 //==============================================================================
 
 std::pair<double, int32_t>
-CSGCell::distance(Position r, Direction u, int32_t on_surface) const
+CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
@@ -781,13 +781,13 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
 DAGCell::DAGCell() : Cell{} {};
 
 std::pair<double, int32_t>
-DAGCell::distance(Position r, Direction u, int32_t on_surface) const
+DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
 {
   // if we've changed direction or we're not on a surface,
   // reset the history and update last direction
-  if (u != simulation::last_dir || on_surface == 0) {
-    simulation::history.reset();
-    simulation::last_dir = u;
+  if (u != p->last_dir_ || on_surface == 0) {
+    p->history_.reset();
+    p->last_dir_ = u;
   }
 
   moab::ErrorCode rval;
@@ -796,7 +796,7 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const
   double dist;
   double pnt[3] = {r.x, r.y, r.z};
   double dir[3] = {u.x, u.y, u.z};
-  rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &simulation::history);
+  rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history_);
   MB_CHK_ERR_CONT(rval);
   int surf_idx;
   if (hit_surf != 0) {

src/dagmc.cpp

@@ -39,13 +39,6 @@ namespace openmc {
 
 const std::string DAGMC_FILENAME = "dagmc.h5m";
 
-namespace simulation {
-
-moab::DagMC::RayHistory history;
-Direction last_dir;
-
-}
-
 namespace model {
 
 moab::DagMC* DAG;

src/geometry.cpp

@@ -378,7 +378,7 @@ BoundaryInfo distance_to_boundary(Particle* p)
     Cell& c {*model::cells[coord.cell]};
 
     // Find the oncoming surface in this cell and the distance to it.
-    auto surface_distance = c.distance(r, u, p->surface_);
+    auto surface_distance = c.distance(r, u, p->surface_, p);
     d_surf = surface_distance.first;
     level_surf_cross = surface_distance.second;
 

src/particle.cpp

@@ -367,7 +367,7 @@ void
 Particle::event_death()
 {
   #ifdef DAGMC
-  if (settings::dagmc) simulation::history.reset();
+  if (settings::dagmc) history_.reset();
   #endif
 
   // Finish particle track output.
@@ -464,7 +464,7 @@ Particle::cross_surface()
     }
 
     Direction u = (surf->bc_ == BC_REFLECT) ?
-      surf->reflect(this->r(), this->u()) :
+      surf->reflect(this->r(), this->u(), this) :
       surf->diffuse_reflect(this->r(), this->u(), this->current_seed());
 
     // Make sure new particle direction is normalized

src/surface.cpp

@@ -184,7 +184,7 @@ Surface::sense(Position r, Direction u) const
 }
 
 Direction
-Surface::reflect(Position r, Direction u) const
+Surface::reflect(Position r, Direction u, Particle* p) const
 {
   // Determine projection of direction onto normal and squared magnitude of
   // normal.
@@ -292,11 +292,11 @@ Direction DAGSurface::normal(Position r) const
   return dir;
 }
 
-Direction DAGSurface::reflect(Position r, Direction u) const
+Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
 {
-  simulation::history.reset_to_last_intersection();
-  simulation::last_dir = Surface::reflect(r, u);
-  return simulation::last_dir;
+  p->history_.reset_to_last_intersection();
+  p->last_dir_ = Surface::reflect(r, u, p);
+  return p->last_dir;
 }
 
 void DAGSurface::to_hdf5(hid_t group_id) const {}