Implement ParticleRay class

include/openmc/particle_data.h

@@ -484,7 +484,7 @@ class GeometryState {
  *   Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
  *   https://doi.org/10.1016/j.anucene.2017.11.032.
  */
-class ParticleData : public GeometryState {
+class ParticleData : virtual public GeometryState {
 private:
   //==========================================================================
   // Data members -- see public: below for descriptions

include/openmc/ray.h

@@ -1,14 +1,14 @@
 #ifndef OPENMC_RAY_H
 #define OPENMC_RAY_H
 
-#include "openmc/particle_data.h"
+#include "openmc/particle.h"
 #include "openmc/position.h"
 
 namespace openmc {
 
 // Base class that implements ray tracing logic, not necessarily through
 // defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
+class Ray : virtual public GeometryState {
 
 public:
   // Initialize from location and direction
@@ -32,6 +32,8 @@ class Ray : public GeometryState {
   // Sets the dist_ variable
   void compute_distance();
 
+  virtual void update_distance();
+
 protected:
   // Records how far the ray has traveled
   double traversal_distance_ {0.0};
@@ -46,5 +48,30 @@ class Ray : public GeometryState {
   unsigned event_counter_ {0};
 };
 
+class ParticleRay : public Ray, public Particle {
+
+public:
+  ParticleRay(
+    Position r, Direction u, ParticleType type_, double time_, double E_)
+    : Ray(r, u)
+  {
+    type() = type_;
+    time() = time_;
+    E() = E_;
+  }
+
+  void on_intersection() override;
+
+  // Sets the dist_ variable
+  void update_distance() override;
+
+  const double& traversal_distance() const { return traversal_distance_; }
+  const double& traversal_mfp() const { return traversal_mfp_; }
+
+protected:
+  // Records how much mean free paths the ray traveled
+  double traversal_mfp_ {0.0};
+};
+
 } // namespace openmc
 #endif // OPENMC_RAY_H

src/ray.cpp

@@ -2,6 +2,8 @@
 
 #include "openmc/error.h"
 #include "openmc/geometry.h"
+#include "openmc/material.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/settings.h"
 
 namespace openmc {
@@ -136,8 +138,8 @@ void Ray::trace()
       cross_lattice(*this, boundary(), settings::verbosity >= 10);
     }
 
-    // Record how far the ray has traveled
-    traversal_distance_ += boundary().distance();
+    update_distance();
+
     inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
 
     // Call the specialized logic for this type of ray. Note that we do not
@@ -165,4 +167,47 @@ void Ray::trace()
   }
 }
 
+void Ray::update_distance()
+{
+  // Record how far the ray has traveled
+  traversal_distance_ += boundary().distance();
+}
+
+void ParticleRay::on_intersection() {}
+
+void ParticleRay::update_distance()
+{
+  Ray::update_distance();
+
+  time() += boundary().distance() / speed();
+
+  // Calculate microscopic and macroscopic cross sections
+  if (material() != MATERIAL_VOID) {
+    if (settings::run_CE) {
+      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
+          density_mult() != density_mult_last()) {
+        // If the material is the same as the last material and the
+        // temperature hasn't changed, we don't need to lookup cross
+        // sections again.
+        model::materials[material()]->calculate_xs(*this);
+      }
+    } else {
+      // Get the MG data; unlike the CE case above, we have to re-calculate
+      // cross sections for every collision since the cross sections may
+      // be angle-dependent
+      data::mg.macro_xs_[material()].calculate_xs(*this);
+
+      // Update the particle's group while we know we are multi-group
+      g_last() = g();
+    }
+  } else {
+    macro_xs().total = 0.0;
+    macro_xs().absorption = 0.0;
+    macro_xs().fission = 0.0;
+    macro_xs().nu_fission = 0.0;
+  }
+
+  traversal_mfp_ += macro_xs().total * boundary().distance();
+}
+
 } // namespace openmc