Fix get_index_in_direction for regular meshes

include/openmc/mesh.h

@@ -196,8 +196,9 @@ class Mesh {
   //! Get bin at a given position in space
   //
   //! \param[in] r Position to get bin for
+  //! \param[in] u Direction of the particle
   //! \return Mesh bin
-  virtual int get_bin(Position r) const = 0;
+  virtual int get_bin(Position r, Direction u) const = 0;
 
   //! Get the number of mesh cells.
   virtual int n_bins() const = 0;
@@ -321,7 +322,7 @@ class StructuredMesh : public Mesh {
 
   virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const;
 
-  int get_bin(Position r) const override;
+  int get_bin(Position r, Direction u) const override;
 
   int n_bins() const override;
 
@@ -360,9 +361,10 @@ class StructuredMesh : public Mesh {
   //! Get mesh indices given a position
   //
   //! \param[in] r Position to get indices for
+  //! \param[in] u Direction of the particle
   //! \param[out] in_mesh Whether position is in mesh
   //! \return Array of mesh indices
-  virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
+  virtual MeshIndex get_indices(Position r, Direction u, bool& in_mesh) const;
 
   //! Get mesh indices corresponding to a mesh bin
   //
@@ -373,8 +375,9 @@ class StructuredMesh : public Mesh {
   //! Get mesh index in a particular direction
   //!
   //! \param[in] r Coordinate to get index for
+  //! \param[in] u Direction of the particle
   //! \param[in] i Direction index
-  virtual int get_index_in_direction(double r, int i) const = 0;
+  virtual int get_index_in_direction(double r, double u, int i) const = 0;
 
   //! Get the coordinate for the mesh grid boundary in the positive direction
   //!
@@ -484,7 +487,7 @@ class RegularMesh : public StructuredMesh {
   RegularMesh(hid_t group);
 
   // Overridden methods
-  int get_index_in_direction(double r, int i) const override;
+  int get_index_in_direction(double r, double u, int i) const override;
 
   virtual std::string get_mesh_type() const override;
 
@@ -537,7 +540,7 @@ class RectilinearMesh : public StructuredMesh {
   RectilinearMesh(hid_t group);
 
   // Overridden methods
-  int get_index_in_direction(double r, int i) const override;
+  int get_index_in_direction(double r, double u, int i) const override;
 
   virtual std::string get_mesh_type() const override;
 
@@ -580,9 +583,10 @@ class CylindricalMesh : public PeriodicStructuredMesh {
   CylindricalMesh(hid_t group);
 
   // Overridden methods
-  virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
+  virtual MeshIndex get_indices(
+    Position r, Direction u, bool& in_mesh) const override;
 
-  int get_index_in_direction(double r, int i) const override;
+  int get_index_in_direction(double r, double u, int i) const override;
 
   virtual std::string get_mesh_type() const override;
 
@@ -645,9 +649,10 @@ class SphericalMesh : public PeriodicStructuredMesh {
   SphericalMesh(hid_t group);
 
   // Overridden methods
-  virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
+  virtual MeshIndex get_indices(
+    Position r, Direction u, bool& in_mesh) const override;
 
-  int get_index_in_direction(double r, int i) const override;
+  int get_index_in_direction(double r, double u, int i) const override;
 
   virtual std::string get_mesh_type() const override;
 
@@ -836,7 +841,7 @@ class MOABMesh : public UnstructuredMesh {
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
 
-  int get_bin(Position r) const override;
+  int get_bin(Position r, Direction u) const override;
 
   int n_bins() const override;
 
@@ -948,8 +953,9 @@ class MOABMesh : public UnstructuredMesh {
   //! Get the mesh cell index for a given position
   //
   //! \param[in] r Position to get index for
+  //! \param[in] u Direction of the particle
   //! \param[in,out] in_mesh Whether position is in the mesh
-  int get_index(const Position& r, bool* in_mesh) const;
+  int get_index(const Position& r, const Direction& u, bool* in_mesh) const;
 
   //! Get the mesh cell index from a bin
   //
@@ -1002,7 +1008,7 @@ class LibMesh : public UnstructuredMesh {
 
   Position sample_element(int32_t bin, uint64_t* seed) const override;
 
-  virtual int get_bin(Position r) const override;
+  virtual int get_bin(Position r, Direction u) const override;
 
   int n_bins() const override;
 
@@ -1092,7 +1098,7 @@ class AdaptiveLibMesh : public LibMesh {
 
   void write(const std::string& filename) const override;
 
-  int get_bin(Position r) const override;
+  int get_bin(Position r, Direction u) const override;
 
 protected:
   // Overridden methods

include/openmc/particle_data.h

@@ -352,6 +352,10 @@ class GeometryState {
   Position& r_born() { return r_born_; }
   const Position& r_born() const { return r_born_; }
 
+  // Direction at birth
+  Direction& u_born() { return u_born_; }
+  const Direction& u_born() const { return u_born_; }
+
   // Coordinates of last collision or reflective/periodic surface
   // crossing for current tallies
   Position& r_last_current() { return r_last_current_; }
@@ -426,6 +430,7 @@ class GeometryState {
   vector<int> cell_last_; //!< coordinates for all levels
 
   Position r_born_;         //!< coordinates at birth
+  Position u_born_;         //!< direction at birth
   Position r_last_current_; //!< coordinates of the last collision or
                             //!< reflective/periodic surface crossing for
                             //!< current tallies

include/openmc/random_ray/flat_source_domain.h

@@ -70,7 +70,7 @@ class FlatSourceDomain {
   }
   int64_t lookup_base_source_region_idx(const GeometryState& p) const;
   SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
-  int64_t lookup_mesh_bin(int64_t sr, Position r) const;
+  int64_t lookup_mesh_bin(int64_t sr, Position r, const Direction& u) const;
   int lookup_mesh_idx(int64_t sr) const;
 
   //----------------------------------------------------------------------------

src/cmfd_solver.cpp

@@ -94,7 +94,7 @@ tensor::Tensor<double> count_bank_sites(
     const auto& site = simulation::source_bank[i];
 
     // determine scoring bin for CMFD mesh
-    int mesh_bin = cmfd::mesh->get_bin(site.r);
+    int mesh_bin = cmfd::mesh->get_bin(site.r, site.u);
 
     // if outside mesh, skip particle
     if (mesh_bin < 0) {

src/eigenvalue.cpp

@@ -660,7 +660,7 @@ void ufs_count_sites()
 double ufs_get_weight(const Particle& p)
 {
   // Determine indices on ufs mesh for current location
-  int mesh_bin = simulation::ufs_mesh->get_bin(p.r());
+  int mesh_bin = simulation::ufs_mesh->get_bin(p.r(), p.u());
   if (mesh_bin < 0) {
     p.write_restart();
     fatal_error("Source site outside UFS mesh!");