Bugfix for short tracks near tally mesh edges

openmc-dev · Feb 28, 2017 · b8ddfac · b8ddfac
1 parent f139ce8
commit b8ddfac
Showing 1 changed file with 26 additions and 10 deletions.
diff --git a/src/tally_filter.F90 b/src/tally_filter.F90
@@ -262,6 +262,9 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
       ! A track can span multiple mesh bins so we need to handle a lot of
       ! intersection logic for tracklength tallies.
 
+      ! ========================================================================
+      ! Determine if the track intersects the tally mesh.
+
       ! Copy the starting and ending coordinates of the particle.  Offset these
       ! just a bit for the purposes of determining if there was an intersection
       ! in case the mesh surfaces coincide with lattice/geometric surfaces which
@@ -296,6 +299,9 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
         end if
       end if
 
+      ! ========================================================================
+      ! Figure out which mesh cell to tally.
+
       ! Copy the un-modified coordinates the particle direction.
       xyz0 = p % last_xyz
       xyz1 = p % coord(1) % xyz
@@ -304,11 +310,20 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
       ! Compute the length of the entire track.
       total_distance = sqrt(sum((xyz1 - xyz0)**2))
 
-      ! If we're looking for the first valid bin, check to see if the particle
-      ! starts inside the mesh.
       if (current_bin == NO_BIN_FOUND) then
+        ! We are looking for the first valid mesh bin.  Check to see if the
+        ! particle starts inside the mesh.
         if (any(ijk0(:m % n_dimension) < 1) &
              .or. any(ijk0(:m % n_dimension) > m % dimension)) then
+          ! The particle does not start in the mesh.  Note that we nudged the
+          ! start and end coordinates by a TINY_BIT each so we will have
+          ! difficulty resolving tracks that are less than 2*TINY_BIT in length.
+          ! If the track is that short, it is also insignificant so we can
+          ! safely ignore it in the tallies.
+          if (total_distance < 2*TINY_BIT) then
+            next_bin = NO_BIN_FOUND
+            return
+          end if
 
           ! The particle does not start in the mesh so keep iterating the ijk0
           ! indices to cross the nearest mesh surface until we've found a valid
@@ -347,14 +362,12 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
           xyz0 = xyz0 + distance * uvw
 
         end if
-      end if
 
-      ! ========================================================================
-      ! If we've already scored some mesh bins, figure out which mesh cell is
-      ! next and where the particle enters that cell.
+      else
+        ! We have already scored some mesh bins for this track.  Pick up where
+        ! we left off and find the next mesh cell that the particle enters.
 
-      if (current_bin /= NO_BIN_FOUND) then
-        ! Get the indices to the last bin.
+        ! Get the indices to the last bin we scored.
         call bin_to_mesh_indices(m, current_bin, ijk0(:m % n_dimension))
 
         ! If the particle track ends in that bin, then we are done.
@@ -401,7 +414,10 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
         end if
       end if
 
-      ! Compute the length of the track segment in this mesh cell.
+      ! ========================================================================
+      ! Compute the length of the track segment in the appropiate mesh cell and
+      ! return.
+
       if (all(ijk0(:m % n_dimension) == ijk1(:m % n_dimension))) then
         ! The track ends in this cell.  Use the particle end location rather
         ! than the mesh surface.
@@ -422,7 +438,7 @@ subroutine get_next_bin_mesh(this, p, estimator, current_bin, next_bin, &
         distance = minval(d(:m % n_dimension))
       end if
 
-      ! Assign the next tally bin and the score
+      ! Assign the next tally bin and the score.
       next_bin = mesh_indices_to_bin(m, ijk0(:m % n_dimension))
       weight = distance / total_distance
     endif