diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 21fc3a472e5..6eadc424299 100644 --- 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