Fix coordinate frame of surface normals in surface tallies and surface sources

include/openmc/geometry.h

@@ -8,6 +8,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/position.h"
 #include "openmc/random_ray/source_region.h" // For hash_combine
 #include "openmc/vector.h"
 
@@ -86,6 +87,23 @@ int check_cell_overlap(GeometryState& p, bool error = true);
 
 int cell_instance_at_level(const GeometryState& p, int level);
 
+//==============================================================================
+//! Rotate a direction from a local coordinate frame into the root frame
+//!
+//! Surface and lattice normals are expressed in the local coordinate frame of
+//! the universe that contains them. Comparing such a normal against the
+//! particle's direction of travel (which is stored in the root frame at
+//! coordinate level zero) requires undoing the rotations that were applied
+//! while descending to \c level.
+//!
+//! \param p A particle whose coordinate levels give the chain of rotations
+//! \param level The level (zero indexed) that \c u is expressed in
+//! \param u A direction in the local frame of \c level
+//! \return The same direction expressed in the root coordinate frame
+//==============================================================================
+
+Direction rotate_to_root(const GeometryState& p, int level, Direction u);
+
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
 //!

include/openmc/particle_data.h

@@ -384,6 +384,14 @@ class GeometryState {
   int& surface() { return surface_; }
   const int& surface() const { return surface_; }
 
+  // Outward unit normal of the surface (or lattice boundary) currently being
+  // crossed, expressed in the root coordinate frame. Set by
+  // score_surface_tally() immediately before the tally filters are evaluated,
+  // so that filters can compare it against the root-frame direction u()
+  // without having to know which coordinate level the surface lives in.
+  Direction& surface_normal() { return surface_normal_; }
+  const Direction& surface_normal() const { return surface_normal_; }
+
   // Surface index based on the current value of the surface_ attribute
   int surface_index() const
   {
@@ -436,6 +444,9 @@ class GeometryState {
   int surface_ {
     SURFACE_NONE}; //!< surface token for surface the particle is currently on
 
+  //! Outward normal of the surface being crossed, in the root coordinate frame
+  Direction surface_normal_ {0.0, 0.0, 1.0};
+
   BoundaryInfo boundary_; //!< Info about the next intersection
 
   int material_ {-1};      //!< index for current material

include/openmc/surface.h

@@ -41,6 +41,21 @@ class Surface {
   unique_ptr<BoundaryCondition> bc_; //!< Boundary condition
   bool surf_source_ {false}; //!< Activate source banking for the surface?
 
+  //! Is this surface used only by cells in the root universe?
+  //!
+  //! evaluate() and normal() work in the local coordinate frame of the
+  //! universe holding the surface. That frame coincides with the root (lab)
+  //! frame only for surfaces in the root universe, so this flag marks the
+  //! surfaces for which a lab-frame position may be passed to them directly.
+  //! Set by finalize_geometry().
+  //!
+  //! Only meaningful for CSG surfaces. finalize_geometry() determines it by
+  //! walking the surfaces named in each cell's region, and DAGCell does not
+  //! report any (it does not override Cell::surfaces()), so a DAGMC surface
+  //! keeps the default here even when its universe is nested below the root
+  //! and transformed. Check geom_type() before relying on this flag.
+  bool root_frame_ {true};
+
   explicit Surface(pugi::xml_node surf_node);
   Surface();
 

src/geometry.cpp

@@ -124,6 +124,23 @@ int cell_instance_at_level(const GeometryState& p, int level)
 
 //==============================================================================
 
+Direction rotate_to_root(const GeometryState& p, int level, Direction u)
+{
+  // Each coordinate level below the root was reached by applying the rotation
+  // matrix of the cell one level above it, so walk back up applying the
+  // inverse of each rotation in turn. Translations are irrelevant here since
+  // they do not affect directions.
+  for (int i = level; i > 0; --i) {
+    if (p.coord(i).rotated()) {
+      const auto& c {*model::cells[p.coord(i - 1).cell()]};
+      u = u.inverse_rotate(c.rotation_);
+    }
+  }
+  return u;
+}
+
+//==============================================================================
+
 bool find_cell_inner(
   GeometryState& p, const NeighborList* neighbor_list, bool verbose)
 {

src/geometry_aux.cpp

@@ -282,6 +282,19 @@ void finalize_geometry()
 
   // Determine number of nested coordinate levels in the geometry
   model::n_coord_levels = maximum_levels(model::root_universe);
+
+  // Flag surfaces that are used outside of the root universe. Their local
+  // coordinate frame does not coincide with the root frame, so a lab-frame
+  // position cannot be handed to Surface::evaluate() or Surface::normal().
+  // Note that DAGCell reports no surfaces, so this leaves DAGMC surfaces at
+  // their default; see the comment on Surface::root_frame_.
+  for (const auto& c : model::cells) {
+    if (c->universe_ == model::root_universe)
+      continue;
+    for (auto token : c->surfaces()) {
+      model::surfaces[std::abs(token) - 1]->root_frame_ = false;
+    }
+  }
 }
 
 //==============================================================================

src/particle.cpp

@@ -344,32 +344,60 @@ void Particle::event_cross_surface()
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
 
+  // The surface or lattice being crossed belongs to the universe at the
+  // coordinate level the boundary search found it on, so its normal is
+  // reported in that level's local frame while the particle direction used to
+  // score surface tallies lives in the root frame. The normal therefore has to
+  // be evaluated at the local position and rotated up into the root frame, and
+  // that has to happen before the crossing is carried out, since crossing
+  // invalidates the coordinate levels.
+  //
+  // Take the level from the boundary rather than from n_coord(). The two are
+  // equal here because of the assignment just above, but reading it from the
+  // boundary keeps this independent of that, and matches where the level came
+  // from originally.
+  int i_surf_level = boundary().coord_level() - 1;
+
   if (boundary().lattice_translation()[0] != 0 ||
       boundary().lattice_translation()[1] != 0 ||
       boundary().lattice_translation()[2] != 0) {
     // Particle crosses lattice boundary
 
-    int i_lattice = coord(boundary().coord_level() - 1).lattice();
-    bool verbose = settings::verbosity >= 10 || trace();
-    cross_lattice(*this, boundary(), verbose);
-    event() = TallyEvent::LATTICE;
+    int i_lattice = coord(i_surf_level).lattice();
 
-    // Score cell to cell partial currents
+    // Determine the lattice boundary normal in the root frame before crossing
+    bool normal_is_valid = false;
+    Direction normal;
     if (!model::active_surface_tallies.empty()) {
       auto& lat {*model::lattices[i_lattice]};
       bool is_valid;
-      Direction normal =
-        lat.get_normal(boundary().lattice_translation(), is_valid);
+      normal = lat.get_normal(boundary().lattice_translation(), is_valid);
       if (is_valid) {
-        normal /= normal.norm();
-        score_surface_tally(*this, model::active_surface_tallies, normal);
+        normal = rotate_to_root(*this, i_surf_level, normal / normal.norm());
+        normal_is_valid = true;
       }
     }
 
+    bool verbose = settings::verbosity >= 10 || trace();
+    cross_lattice(*this, boundary(), verbose);
+    event() = TallyEvent::LATTICE;
+
+    // Score cell to cell partial currents
+    if (normal_is_valid) {
+      score_surface_tally(*this, model::active_surface_tallies, normal);
+    }
+
   } else {
 
     const auto& surf {*model::surfaces[surface_index()].get()};
 
+    // Determine the surface normal in the root frame before crossing
+    Direction normal;
+    if (!model::active_surface_tallies.empty()) {
+      normal = surf.normal(coord(i_surf_level).r());
+      normal = rotate_to_root(*this, i_surf_level, normal / normal.norm());
+    }
+
     // Particle crosses surface
     // If BC, add particle to surface source before crossing surface
     if (surf.surf_source_ && surf.bc_) {
@@ -387,8 +415,6 @@ void Particle::event_cross_surface()
 
     // 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);
     }
   }

src/plot.cpp

@@ -1825,12 +1825,7 @@ void PhongRay::on_intersection()
 
     // Need to apply rotations to find the normal vector in
     // the base level universe's coordinate system.
-    for (int lev = surf_level - 1; lev >= 0; --lev) {
-      if (coord(lev + 1).rotated()) {
-        const Cell& c {*model::cells[coord(lev).cell()]};
-        normal = normal.inverse_rotate(c.rotation_);
-      }
-    }
+    normal = rotate_to_root(*this, surf_level, normal);
 
     // use the normal opposed to the ray direction
     if (normal.dot(u()) > 0.0) {

src/source.cpp

@@ -554,11 +554,19 @@ SourceSite FileSource::sample(uint64_t* seed) const
   // surface containing the source site, determine the signed half-space from
   // the particle direction. Otherwise, ignore the surface ID and allow the
   // normal cell search to locate the particle.
+  //
+  // The site position and direction are in the root coordinate frame, whereas
+  // evaluate() and normal() work in the local frame of the universe holding
+  // the surface. The half-space can therefore only be recovered here for
+  // surfaces in the root universe; for any other surface the frames differ by
+  // an unknown transform (a universe may be filled in several places with
+  // different rotations, so the transform cannot be recovered from the site
+  // alone) and the surface ID is dropped rather than signed incorrectly.
   if (site.surf_id != SURFACE_NONE) {
     auto it = model::surface_map.find(std::abs(site.surf_id));
     if (it != model::surface_map.end()) {
       const auto& surf = *model::surfaces[it->second];
-      if (surf.geom_type() == GeometryType::CSG &&
+      if (surf.geom_type() == GeometryType::CSG && surf.root_frame_ &&
           std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) {
         int surf_id = std::abs(site.surf_id);
         site.surf_id =

src/tallies/filter_musurface.cpp

@@ -3,18 +3,18 @@
 #include <cmath> // for abs, copysign
 
 #include "openmc/search.h"
-#include "openmc/surface.h"
 #include "openmc/tallies/tally_scoring.h"
 
 namespace openmc {
 
 void MuSurfaceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  // Get surface normal (and make sure it is a unit vector)
-  const auto surf {model::surfaces[p.surface_index()].get()};
-  auto n = surf->normal(p.r());
-  n /= n.norm();
+  // Use the normal recorded for the crossing being scored. It is already a
+  // unit vector expressed in the root coordinate frame, which is the frame
+  // p.u() is in -- recomputing it here from the surface would give the normal
+  // in the local frame of whichever universe holds the surface.
+  Direction n = p.surface_normal();
 
   // Determine whether normal should be pointing in or out
   if (p.surface() < 0)

src/tallies/tally_scoring.cpp

@@ -2662,6 +2662,11 @@ void score_surface_tally(
 {
   double wgt = p.wgt_last();
 
+  // Make the normal available to filters (e.g. MuSurfaceFilter) that need it.
+  // The caller is responsible for supplying it in the root coordinate frame so
+  // that it can be compared directly against p.u().
+  p.surface_normal() = normal;
+
   double mu = std::clamp(p.u().dot(normal), -1.0, 1.0);
 
   // Sign for net current: +1 if crossing outward (in direction of normal),

tests/unit_tests/test_filter_musurface.py

@@ -1,3 +1,5 @@
+import math
+
 import openmc
 
 
@@ -36,3 +38,49 @@ def test_musurface(run_in_tmpdir):
         assert element == 0.0
 
 
+def test_musurface_rotated_universe(run_in_tmpdir):
+    """MuSurfaceFilter uses the surface normal in the root coordinate frame.
+
+    The plane lives in a universe filled into a cell rotated 45 degrees about
+    z, so its normal in the plane's own frame is (1, 0, 0) while the particle
+    direction is stored in the root frame. Binning must use the root-frame
+    normal, giving mu = cos(45 deg) rather than 1.
+    """
+    openmc.reset_auto_ids()
+
+    xplane = openmc.XPlane(0.0)
+    inner1 = openmc.Cell(region=-xplane)
+    inner2 = openmc.Cell(region=+xplane)
+    inner_univ = openmc.Universe(cells=[inner1, inner2])
+
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    root_cell = openmc.Cell(region=-sph, fill=inner_univ)
+    root_cell.rotation = (0.0, 0.0, 45.0)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root_cell])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    # 20 equal-width bins from -1 to 1; cos(45 deg) = 0.7071 falls in [0.7, 0.8)
+    tally = openmc.Tally()
+    tally.filters = [
+        openmc.MuSurfaceFilter(20),
+        openmc.SurfaceFilter([xplane]),
+    ]
+    tally.scores = ['current']
+    model.tallies = [tally]
+
+    model.run(apply_tally_results=True)
+    current_mu = tally.mean.ravel()
+
+    expected_bin = int((math.cos(math.radians(45.0)) + 1.0) / 0.1)
+    assert current_mu[expected_bin] == 1.0
+    assert current_mu.sum() == 1.0