From 6cd39073b3523768748138f90be9e31d19cfdb52 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Mon, 23 Mar 2026 17:16:59 +0200 Subject: [PATCH] Fix surface tally when crossing lattice (#3895) --- include/openmc/lattice.h | 14 ++++ include/openmc/tallies/tally_scoring.h | 4 +- src/lattice.cpp | 105 +++++++++++++++++++++++++ src/particle.cpp | 34 ++++++-- src/tallies/tally_scoring.cpp | 10 +-- 5 files changed, 153 insertions(+), 14 deletions(-) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index f87d28b21..ca40bbc2a 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -113,6 +113,14 @@ public: virtual Position get_local_position( Position r, const array& 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& 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 @@ public: Position get_local_position( Position r, const array& i_xyz) const override; + Direction get_normal( + const array& i_xyz, bool& is_valid) const override; + int32_t& offset(int map, const array& i_xyz) override; int32_t offset(int map, int indx) const override; @@ -268,6 +279,9 @@ public: Position get_local_position( Position r, const array& i_xyz) const override; + Direction get_normal( + const array& i_xyz, bool& is_valid) const override; + bool is_valid_index(int indx) const override; int32_t& offset(int map, const array& i_xyz) override; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 29b3ec6e5..d1aed2831 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -111,9 +111,9 @@ void score_meshsurface_tally(Particle& p, const vector& 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& tallies, const Surface& surf); + Particle& p, const vector& tallies, const Direction& normal); //! Score the pulse-height tally //! This is triggered at the end of every particle history diff --git a/src/lattice.cpp b/src/lattice.cpp index 92d451f61..e799a340e 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -340,6 +340,26 @@ Position RectLattice::get_local_position( //============================================================================== +Direction RectLattice::get_normal( + const array& 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& 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& 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}; diff --git a/src/particle.cpp b/src/particle.cpp index 0a0635598..11747e0cc 100644 --- a/src/particle.cpp +++ b/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()) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 4210b034a..d5af99382 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2656,19 +2656,19 @@ void score_meshsurface_tally(Particle& p, const vector& tallies) } void score_surface_tally( - Particle& p, const vector& tallies, const Surface& surf) + Particle& p, const vector& 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;