From 521964149f39d66643d40a546eb254a52f4403c1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 28 Oct 2019 11:15:53 -0400 Subject: [PATCH] Fix tally mesh bug for very short tracks --- src/mesh.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index ac02b679e..d6da64bb7 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -523,6 +523,14 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } r1 = r; + // The TINY_BIT offsets above mean that the preceding logic cannot always find + // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just + // assume the track lies in only one mesh bin. These tracks are very short so + // any error caused by this assumption will be small. + if (total_distance < 2*TINY_BIT) { + for (int i = 0; i < n; ++i) ijk0[i] = ijk1[i]; + } + // ======================================================================== // Find which mesh cells are traversed and the length of each traversal. @@ -898,6 +906,14 @@ void RectilinearMesh::bins_crossed(const Particle* p, std::vector& bins, } r1 = r; + // The TINY_BIT offsets above mean that the preceding logic cannot always find + // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just + // assume the track lies in only one mesh bin. These tracks are very short so + // any error caused by this assumption will be small. + if (total_distance < 2*TINY_BIT) { + for (int i = 0; i < 3; ++i) ijk0[i] = ijk1[i]; + } + // ======================================================================== // Find which mesh cells are traversed and the length of each traversal.