Merge pull request #809 from smharper/mesh_error

Bugfix for short tracks near tally mesh edges
This commit is contained in:
Paul Romano 2017-03-01 06:21:42 -06:00 committed by GitHub
commit e22bb71ad7

View file

@ -262,6 +262,9 @@ contains
! 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 @@ contains
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 @@ contains
! 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 @@ contains
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 @@ contains
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 @@ contains
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