Fix surface tally when crossing lattice (#3895)

Author: GuySten

include/openmc/lattice.h

@@ -113,6 +113,14 @@ class Lattice {
   virtual Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const = 0;
 
+  //! \brief get the normal of the lattice surface crossing
+  //! \param[in] i_xyz The indices for the lattice translation.
+  //! \param[out] is_valid is the lattice translation correspond to a valid
+  //! surface. \return The surface normal corresponding to the lattice
+  //! translation.
+  virtual Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const = 0;
+
   //! \brief Check flattened lattice index.
   //! \param indx The index for a lattice tile.
   //! \return true if the given index fit within the lattice bounds.  False
@@ -223,6 +231,9 @@ class RectLattice : public Lattice {
   Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const override;
 
+  Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const override;
+
   int32_t& offset(int map, const array<int, 3>& i_xyz) override;
 
   int32_t offset(int map, int indx) const override;
@@ -268,6 +279,9 @@ class HexLattice : public Lattice {
   Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const override;
 
+  Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const override;
+
   bool is_valid_index(int indx) const override;
 
   int32_t& offset(int map, const array<int, 3>& i_xyz) override;

include/openmc/tallies/tally_scoring.h

@@ -111,9 +111,9 @@ void score_meshsurface_tally(Particle& p, const vector<int>& tallies);
 //
 //! \param p The particle being tracked
 //! \param tallies A vector of the indices of the tallies to score to
-//! \param surf The surface being crossed
+//! \param normal The normal of the surface being crossed
 void score_surface_tally(
-  Particle& p, const vector<int>& tallies, const Surface& surf);
+  Particle& p, const vector<int>& tallies, const Direction& normal);
 
 //! Score the pulse-height tally
 //! This is triggered at the end of every particle history

src/lattice.cpp

@@ -340,6 +340,26 @@ Position RectLattice::get_local_position(
 
 //==============================================================================
 
+Direction RectLattice::get_normal(
+  const array<int, 3>& i_xyz, bool& is_valid) const
+{
+  is_valid = false;
+  Direction dir = {0.0, 0.0, 0.0};
+  if ((std::abs(i_xyz[0]) == 1) && (i_xyz[1] == 0) && (i_xyz[2] == 0)) {
+    is_valid = true;
+    dir[0] = std::copysign(1.0, i_xyz[0]);
+  } else if ((i_xyz[0] == 0) && (std::abs(i_xyz[1]) == 1) && (i_xyz[2] == 0)) {
+    is_valid = true;
+    dir[1] = std::copysign(1.0, i_xyz[1]);
+  } else if ((i_xyz[0] == 0) && (i_xyz[1] == 0) && (std::abs(i_xyz[2]) == 1)) {
+    is_valid = true;
+    dir[2] = std::copysign(1.0, i_xyz[2]);
+  }
+  return dir;
+}
+
+//==============================================================================
+
 int32_t& RectLattice::offset(int map, const array<int, 3>& i_xyz)
 {
   return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map +
@@ -986,6 +1006,91 @@ Position HexLattice::get_local_position(
 
 //==============================================================================
 
+Direction HexLattice::get_normal(
+  const array<int, 3>& i_xyz, bool& is_valid) const
+{
+  // Short description of the direction vectors used here.  The beta, gamma, and
+  // delta vectors point towards the flat sides of each hexagonal tile.
+  // Y - orientation:
+  //   basis0 = (1, 0)
+  //   basis1 = (-1/sqrt(3), 1)   = +120 degrees from basis0
+  //   beta   = (sqrt(3)/2, 1/2)  = +30 degrees from basis0
+  //   gamma  = (sqrt(3)/2, -1/2) = -60 degrees from beta
+  //   delta  = (0, 1)            = +60 degrees from beta
+  // X - orientation:
+  //   basis0 = (1/sqrt(3), -1)
+  //   basis1 = (0, 1)            = +120 degrees from basis0
+  //   beta   = (1, 0)            = +30 degrees from basis0
+  //   gamma  = (1/2, -sqrt(3)/2) = -60 degrees from beta
+  //   delta  = (1/2, sqrt(3)/2)  = +60 degrees from beta
+
+  is_valid = false;
+  Direction dir = {0.0, 0.0, 0.0};
+  if ((i_xyz[0] == 0) && (i_xyz[1] == 0) && (std::abs(i_xyz[2]) == 1)) {
+    is_valid = true;
+    dir[2] = std::copysign(1.0, i_xyz[2]);
+  } else if ((i_xyz[2] == 0) &&
+             std::max({std::abs(i_xyz[0]), std::abs(i_xyz[1]),
+               std::abs(i_xyz[0] + i_xyz[1])}) == 1) {
+    is_valid = true;
+    // beta direction
+    if ((i_xyz[0] == 1) && (i_xyz[1] == 0)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.5 * std::sqrt(3.0);
+        dir[1] = 0.5;
+      } else {
+        dir[0] = 1.0;
+        dir[1] = 0.0;
+      }
+    } else if ((i_xyz[0] == -1) && (i_xyz[1] == 0)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = -0.5 * std::sqrt(3.0);
+        dir[1] = -0.5;
+      } else {
+        dir[0] = -1.0;
+        dir[1] = 0.0;
+      }
+      // gamma direction
+    } else if ((i_xyz[0] == 1) && (i_xyz[1] == -1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.5 * std::sqrt(3.0);
+        dir[1] = -0.5;
+      } else {
+        dir[0] = 0.5;
+        dir[1] = -0.5 * std::sqrt(3.0);
+      }
+    } else if ((i_xyz[0] == -1) && (i_xyz[1] == 1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = -0.5 * std::sqrt(3.0);
+        dir[1] = 0.5;
+      } else {
+        dir[0] = -0.5;
+        dir[1] = 0.5 * std::sqrt(3.0);
+      }
+      // delta direction
+    } else if ((i_xyz[0] == 0) && (i_xyz[1] == 1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.0;
+        dir[1] = 1.0;
+      } else {
+        dir[0] = 0.5;
+        dir[1] = 0.5 * std::sqrt(3.0);
+      }
+    } else if ((i_xyz[0] == 0) && (i_xyz[1] == -1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.0;
+        dir[1] = -1.0;
+      } else {
+        dir[0] = -0.5;
+        dir[1] = -0.5 * std::sqrt(3.0);
+      }
+    }
+  }
+  return dir;
+}
+
+//==============================================================================
+
 bool HexLattice::is_valid_index(int indx) const
 {
   int nx {2 * n_rings_ - 1};

src/particle.cpp

@@ -14,6 +14,7 @@
 #include "openmc/error.h"
 #include "openmc/geometry.h"
 #include "openmc/hdf5_interface.h"
+#include "openmc/lattice.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
 #include "openmc/mgxs_interface.h"
@@ -302,8 +303,6 @@ void Particle::event_cross_surface()
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
 
-  const auto& surf {*model::surfaces[surface_index()].get()};
-
   if (boundary().lattice_translation()[0] != 0 ||
       boundary().lattice_translation()[1] != 0 ||
       boundary().lattice_translation()[2] != 0) {
@@ -312,7 +311,23 @@ void Particle::event_cross_surface()
     bool verbose = settings::verbosity >= 10 || trace();
     cross_lattice(*this, boundary(), verbose);
     event() = TallyEvent::LATTICE;
+
+    // Score cell to cell partial currents
+    if (!model::active_surface_tallies.empty()) {
+      auto& lat {*model::lattices[lowest_coord().lattice()]};
+      bool is_valid;
+      Direction normal =
+        lat.get_normal(boundary().lattice_translation(), is_valid);
+      if (is_valid) {
+        normal /= normal.norm();
+        score_surface_tally(*this, model::active_surface_tallies, normal);
+      }
+    }
+
   } else {
+
+    const auto& surf {*model::surfaces[surface_index()].get()};
+
     // Particle crosses surface
     // If BC, add particle to surface source before crossing surface
     if (surf.surf_source_ && surf.bc_) {
@@ -327,10 +342,13 @@ void Particle::event_cross_surface()
       apply_weight_windows(*this);
     }
     event() = TallyEvent::SURFACE;
-  }
-  // Score cell to cell partial currents
-  if (!model::active_surface_tallies.empty()) {
-    score_surface_tally(*this, model::active_surface_tallies, surf);
+
+    // Score cell to cell partial currents
+    if (!model::active_surface_tallies.empty()) {
+      Direction normal = surf.normal(r());
+      normal /= normal.norm();
+      score_surface_tally(*this, model::active_surface_tallies, normal);
+    }
   }
 }
 
@@ -681,7 +699,9 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   // with a mesh boundary
 
   if (!model::active_surface_tallies.empty()) {
-    score_surface_tally(*this, model::active_surface_tallies, surf);
+    Direction normal = surf.normal(r());
+    normal /= normal.norm();
+    score_surface_tally(*this, model::active_surface_tallies, normal);
   }
 
   if (!model::active_meshsurf_tallies.empty()) {

src/tallies/tally_scoring.cpp

@@ -2656,19 +2656,19 @@ void score_meshsurface_tally(Particle& p, const vector<int>& tallies)
 }
 
 void score_surface_tally(
-  Particle& p, const vector<int>& tallies, const Surface& surf)
+  Particle& p, const vector<int>& tallies, const Direction& normal)
 {
   double wgt = p.wgt_last();
 
+  double mu = std::clamp(p.u().dot(normal), -1.0, 1.0);
+
   // Sign for net current: +1 if crossing outward (in direction of normal),
   // -1 if crossing inward
-  double current_sign = (p.surface() > 0) ? 1.0 : -1.0;
+  double current_sign = std::copysign(1.0, mu);
 
   // Determine absolute cosine of angle between particle direction and surface
   // normal, needed for the surface-crossing flux estimator.
-  auto n = surf.normal(p.r());
-  n /= n.norm();
-  double abs_mu = std::min(std::abs(p.u().dot(n)), 1.0);
+  double abs_mu = std::abs(mu);
   if (abs_mu < settings::surface_grazing_cutoff)
     abs_mu = settings::surface_grazing_ratio * settings::surface_grazing_cutoff;