Merge pull request #1255 from pshriwise/dagmc_rayhist

Adding RayHistory to DagMC

include/openmc/cell.h

@@ -10,9 +10,7 @@
 
 #include "hdf5.h"
 #include "pugixml.hpp"
-#ifdef DAGMC
-#include "DagMC.hpp"
-#endif
+#include "dagmc.h"
 
 #include "openmc/constants.h"
 #include "openmc/neighbor_list.h"

include/openmc/dagmc.h

@@ -9,11 +9,19 @@ extern "C" const bool dagmc_enabled;
 #ifdef DAGMC
 
 #include "DagMC.hpp"
-#include "openmc/cell.h"
-#include "openmc/surface.h"
+#include "openmc/xml_interface.h"
+#include "openmc/position.h"
 
 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/surface.h

@@ -12,10 +12,7 @@
 
 #include "openmc/constants.h"
 #include "openmc/position.h"
-
-#ifdef DAGMC
-#include "DagMC.hpp"
-#endif
+#include "dagmc.h"
 
 namespace openmc {
 
@@ -82,7 +79,7 @@ public:
   //! \param[in] r The point at which the ray is incident.
   //! \param[in] u Incident direction of the ray
   //! \return Outgoing direction of the ray
-  Direction reflect(Position r, Direction u) const;
+  virtual Direction reflect(Position r, Direction u) const;
 
   //! Evaluate the equation describing the surface.
   //!
@@ -136,6 +133,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;
   //! Get the bounding box of this surface.
   BoundingBox bounding_box() const;
 

src/cell.cpp

@@ -1,3 +1,4 @@
+
 #include "openmc/cell.h"
 
 #include <cmath>
@@ -8,6 +9,7 @@
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
+#include "openmc/dagmc.h"
 #include "openmc/error.h"
 #include "openmc/geometry.h"
 #include "openmc/hdf5_interface.h"
@@ -614,19 +616,28 @@ DAGCell::DAGCell() : Cell{} {};
 std::pair<double, int32_t>
 DAGCell::distance(Position r, Direction u, int32_t on_surface) 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;
+  }
+
   moab::ErrorCode rval;
   moab::EntityHandle vol = dagmc_ptr_->entity_by_index(3, dag_index_);
   moab::EntityHandle hit_surf;
   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);
+  rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &simulation::history);
   MB_CHK_ERR_CONT(rval);
   int surf_idx;
   if (hit_surf != 0) {
     surf_idx = dagmc_ptr_->index_by_handle(hit_surf);
-  } else {  // indicate that particle is lost
+  } else {
+    // indicate that particle is lost
     surf_idx = -1;
+    dist = INFINITY;
   }
 
   return {dist, surf_idx};

src/dagmc.cpp

@@ -9,6 +9,7 @@
 #include "openmc/material.h"
 #include "openmc/string_utils.h"
 #include "openmc/settings.h"
+#include "openmc/surface.h"
 
 #ifdef DAGMC
 
@@ -38,6 +39,14 @@ const std::string DAGMC_FILENAME = "dagmc.h5m";
 
 namespace openmc {
 
+
+namespace simulation {
+
+moab::DagMC::RayHistory history;
+Direction last_dir;
+
+}
+
 namespace model {
 
 moab::DagMC* DAG;

src/particle.cpp

@@ -372,6 +372,10 @@ Particle::transport()
     }
   }
 
+  #ifdef DAGMC
+  if (settings::dagmc) simulation::history.reset();
+  #endif
+
   // Finish particle track output.
   if (write_track_) {
     write_particle_track(*this);
@@ -460,12 +464,14 @@ Particle::cross_surface()
     // If a reflective surface is coincident with a lattice or universe
     // boundary, it is necessary to redetermine the particle's coordinates in
     // the lower universes.
-
-    n_coord_ = 1;
-    if (!find_cell(this, true)) {
-      this->mark_as_lost("Couldn't find particle after reflecting from surface "
-        + std::to_string(surf->id_) + ".");
-      return;
+    // (unless we're using a dagmc model, which has exactly one universe)
+    if (!settings::dagmc) {
+      n_coord_ = 1;
+      if (!find_cell(this, true)) {
+        this->mark_as_lost("Couldn't find particle after reflecting from surface "
+                           + std::to_string(surf->id_) + ".");
+        return;
+      }
     }
 
     // Set previous coordinate going slightly past surface crossing

src/surface.cpp

@@ -6,6 +6,7 @@
 #include <utility>
 
 #include "openmc/error.h"
+#include "openmc/dagmc.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/settings.h"
 #include "openmc/string_utils.h"
@@ -255,13 +256,19 @@ DAGSurface::distance(Position r, Direction u, bool coincident) const
 Direction DAGSurface::normal(Position r) const
 {
   moab::ErrorCode rval;
-  Direction u;
   moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
   double pnt[3] = {r.x, r.y, r.z};
-  double dir[3] = {u.x, u.y, u.z};
+  double dir[3];
   rval = dagmc_ptr_->get_angle(surf, pnt, dir);
   MB_CHK_ERR_CONT(rval);
-  return u;
+  return dir;
+}
+
+Direction DAGSurface::reflect(Position r, Direction u) const
+{
+  simulation::history.reset_to_last_intersection();
+  simulation::last_dir = Surface::reflect(r, u);
+  return simulation::last_dir;
 }
 
 BoundingBox DAGSurface::bounding_box() const

tests/regression_tests/dagmc/refl/results_true.dat

@@ -0,0 +1,5 @@
+k-combined:
+2.053871E+00 7.718195E-03
+tally 1:
+1.171003E+01
+2.864565E+01

tests/regression_tests/dagmc/refl/test.py

@@ -0,0 +1,40 @@
+import openmc
+import openmc.capi
+from openmc.stats import Box
+
+import pytest
+from tests.testing_harness import PyAPITestHarness
+
+pytestmark = pytest.mark.skipif(
+    not openmc.capi._dagmc_enabled(),
+    reason="DAGMC CAD geometry is not enabled.")
+
+class UWUWTest(PyAPITestHarness):
+
+    def _build_inputs(self):
+        model = openmc.model.Model()
+
+        # settings
+        model.settings.batches = 5
+        model.settings.inactive = 0
+        model.settings.particles = 100
+
+        source = openmc.Source(space=Box([-4, -4, -4],
+                                         [ 4,  4,  4]))
+        model.settings.source = source
+
+        model.settings.dagmc = True
+
+        model.settings.export_to_xml()
+
+        # tally
+        tally = openmc.Tally()
+        tally.scores = ['total']
+        tally.filters = [openmc.CellFilter(1)]
+        model.tallies = [tally]
+
+        model.tallies.export_to_xml()
+
+def test_refl():
+    harness = UWUWTest('statepoint.5.h5')
+    harness.main()

tests/regression_tests/dagmc/uwuw/inputs_true.dat

@@ -0,0 +1,23 @@
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>100</particles>
+  <batches>5</batches>
+  <inactive>0</inactive>
+  <source strength="1.0">
+    <space type="box">
+      <parameters>-4 -4 -4 4 4 4</parameters>
+    </space>
+  </source>
+  <dagmc>true</dagmc>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>1</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <scores>total</scores>
+  </tally>
+</tallies>

tests/regression_tests/dagmc/uwuw/results_true.dat

@@ -0,0 +1 @@
+../legacy/results_true.dat

tests/regression_tests/dagmc/uwuw/test.py

@@ -1,7 +1,6 @@
 import openmc
 import openmc.capi
 from openmc.stats import Box
-from openmc.material import Materials
 
 import pytest
 from tests.testing_harness import PyAPITestHarness

tests/regression_tests/uwuw/results_true.dat

@@ -1 +0,0 @@
-../dagmc/results_true.dat