From 47ee328ca17c9080e723b001b487e9986a73244e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 1 Nov 2011 10:51:01 -0400 Subject: [PATCH] Surface current tallies are now called independent of other tallies. --- src/physics.f90 | 3 +- src/tally.f90 | 370 ++++++++++++++++++++++++------------------------ 2 files changed, 190 insertions(+), 183 deletions(-) diff --git a/src/physics.f90 b/src/physics.f90 index 840c16c623..fa693702b3 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -13,7 +13,7 @@ module physics use mcnp_random, only: rang use output, only: write_message, print_particle use particle_header, only: Particle - use tally, only: score_tally + use tally, only: score_tally, score_surface_current use search, only: binary_search use string, only: int_to_str @@ -412,6 +412,7 @@ contains if (tallies_on) then call score_tally(p, scattered, fissioned) + call score_surface_current(p) end if ! Reset number of particles banked during collision diff --git a/src/tally.f90 b/src/tally.f90 index be09cd69f3..10c1dbba06 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -330,11 +330,8 @@ contains do i = 1, n_tallies t => tallies(i) - ! Handle surface current tallies separately - if (t % surface_current) then - call score_surface_current(p, t) - cycle - end if + ! Surface current tallies are treated separately + if (t % surface_current) cycle ! ======================================================================= ! DETERMINE SCORING BIN COMBINATION @@ -573,16 +570,17 @@ contains end subroutine score_tally !=============================================================================== -! SCORE_SURFACE_CURRENT +! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually +! determining which mesh surfaces were crossed !=============================================================================== - subroutine score_surface_current(p, t) + subroutine score_surface_current(p) type(Particle), pointer :: p - type(TallyObject), pointer :: t integer :: i ! loop indices integer :: j ! loop indices + integer :: k ! loop indices integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates integer :: n_cross ! number of surface crossings @@ -598,203 +596,211 @@ contains logical :: x_same ! same starting/ending x index (i) logical :: y_same ! same starting/ending y index (j) logical :: z_same ! same starting/ending z index (k) + type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() - ! Get pointer to mesh - m => meshes(t % mesh) + do i = 1, n_tallies + ! Copy starting and ending location of particle + xyz0 = p % last_xyz + xyz1 = p % xyz - ! Copy starting and ending location of particle - xyz0 = p % last_xyz - xyz1 = p % xyz + ! Get pointer to tally and mesh + t => tallies(i) + m => meshes(t % mesh) - ! Determine indices for starting and ending location - call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) - call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) + ! Skip non-surface-current tallies + if (.not. t % surface_current) cycle - ! Check to make sure start or end is in mesh - if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) return + ! Determine indices for starting and ending location + call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) + call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) - ! Calculate number of surface crossings - n_cross = sum(abs(ijk1 - ijk0)) - if (n_cross == 0) return + ! Check to make sure start or end is in mesh + if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) cycle - ! Copy particle's direction - uvw = p % uvw + ! Calculate number of surface crossings + n_cross = sum(abs(ijk1 - ijk0)) + if (n_cross == 0) cycle - ! ========================================================================== - ! SPECIAL CASES WHERE TWO INDICES ARE THE SAME + ! Copy particle's direction + uvw = p % uvw - x_same = (ijk0(1) == ijk1(1)) - y_same = (ijk0(2) == ijk1(2)) - z_same = (ijk0(3) == ijk1(3)) + ! ======================================================================= + ! SPECIAL CASES WHERE TWO INDICES ARE THE SAME - if (x_same .and. y_same) then - ! Only z crossings - if (uvw(3) > 0) then - do i = ijk0(3), ijk1(3) - 1 - ijk0(3) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - else - do i = ijk0(3) - 1, ijk1(3), -1 - ijk0(3) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_TOP - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - end if - return - elseif (x_same .and. z_same) then - ! Only y crossings - if (uvw(2) > 0) then - do i = ijk0(2), ijk1(2) - 1 - ijk0(2) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - else - do i = ijk0(2) - 1, ijk1(2), -1 - ijk0(2) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - end if - return - elseif (y_same .and. z_same) then - ! Only x crossings - if (uvw(1) > 0) then - do i = ijk0(1), ijk1(1) - 1 - ijk0(1) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - else - do i = ijk0(1) - 1, ijk1(1), -1 - ijk0(1) = i - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - end do - end if - return - end if + x_same = (ijk0(1) == ijk1(1)) + y_same = (ijk0(2) == ijk1(2)) + z_same = (ijk0(3) == ijk1(3)) - ! ========================================================================== - ! GENERIC CASE - - ! Bounding coordinates - do i = 1, 3 - if (uvw(i) > 0) then - xyz_cross(i) = m % origin(i) + ijk0(i) * m % width(i) - else - xyz_cross(i) = m % origin(i) + (ijk0(i) - 1) * m % width(i) - end if - end do - - do i = 1, n_cross - ! Reset scoring bin index - score_index = 0 - - ! Calculate distance to each bounding surface. We need to treat special - ! case where the cosine of the angle is zero since this would result in a - ! divide-by-zero. - - do j = 1, 3 - if (uvw(j) == 0) then - d(j) = INFINITY + if (x_same .and. y_same) then + ! Only z crossings + if (uvw(3) > 0) then + do j = ijk0(3), ijk1(3) - 1 + ijk0(3) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do else - d(j) = (xyz_cross(j) - xyz0(j))/uvw(j) + do j = ijk0(3) - 1, ijk1(3), -1 + ijk0(3) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_TOP + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do + end if + cycle + elseif (x_same .and. z_same) then + ! Only y crossings + if (uvw(2) > 0) then + do j = ijk0(2), ijk1(2) - 1 + ijk0(2) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do + else + do j = ijk0(2) - 1, ijk1(2), -1 + ijk0(2) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do + end if + cycle + elseif (y_same .and. z_same) then + ! Only x crossings + if (uvw(1) > 0) then + do j = ijk0(1), ijk1(1) - 1 + ijk0(1) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do + else + do j = ijk0(1) - 1, ijk1(1), -1 + ijk0(1) = j + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT + call add_to_score(t % scores(score_index, 1), p % last_wgt) + end if + end do + end if + cycle + end if + + ! ======================================================================= + ! GENERIC CASE + + ! Bounding coordinates + do j = 1, 3 + if (uvw(j) > 0) then + xyz_cross(j) = m % origin(j) + ijk0(j) * m % width(j) + else + xyz_cross(j) = m % origin(j) + (ijk0(j) - 1) * m % width(j) end if end do - ! Determine the closest bounding surface of the mesh cell by calculating - ! the minimum distance + do k = 1, n_cross + ! Reset scoring bin index + score_index = 0 - distance = minval(d) - - ! Now use the minimum distance and diretion of the particle to determine - ! which surface was crossed - - if (distance == d(1)) then - if (uvw(1) > 0) then - ! Crossing into right mesh cell -- this is treated as outgoing - ! current from (i,j,k) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT + ! Calculate distance to each bounding surface. We need to treat + ! special case where the cosine of the angle is zero since this would + ! result in a divide-by-zero. + + do j = 1, 3 + if (uvw(j) == 0) then + d(j) = INFINITY + else + d(j) = (xyz_cross(j) - xyz0(j))/uvw(j) end if - ijk0(1) = ijk0(1) + 1 - xyz_cross(1) = xyz_cross(1) + m % width(1) - else - ! Crossing into left mesh cell -- this is treated as incoming - ! current in (i-1,j,k) - ijk0(1) = ijk0(1) - 1 - xyz_cross(1) = xyz_cross(1) - m % width(1) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT + end do + + ! Determine the closest bounding surface of the mesh cell by + ! calculating the minimum distance + + distance = minval(d) + + ! Now use the minimum distance and diretion of the particle to + ! determine which surface was crossed + + if (distance == d(1)) then + if (uvw(1) > 0) then + ! Crossing into right mesh cell -- this is treated as outgoing + ! current from (i,j,k) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT + end if + ijk0(1) = ijk0(1) + 1 + xyz_cross(1) = xyz_cross(1) + m % width(1) + else + ! Crossing into left mesh cell -- this is treated as incoming + ! current in (i-1,j,k) + ijk0(1) = ijk0(1) - 1 + xyz_cross(1) = xyz_cross(1) - m % width(1) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT + end if + end if + elseif (distance == d(2)) then + if (uvw(2) > 0) then + ! Crossing into front mesh cell -- this is treated as outgoing + ! current in (i,j,k) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT + end if + ijk0(2) = ijk0(2) + 1 + xyz_cross(2) = xyz_cross(2) + m % width(2) + else + ! Crossing into back mesh cell -- this is treated as incoming + ! current in (i,j-1,k) + ijk0(2) = ijk0(2) - 1 + xyz_cross(2) = xyz_cross(2) - m % width(2) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT + end if + end if + else if (distance == d(3)) then + if (uvw(3) > 0) then + ! Crossing into top mesh cell -- this is treated as outgoing + ! current in (i,j,k) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP + end if + ijk0(3) = ijk0(3) + 1 + xyz_cross(3) = xyz_cross(3) + m % width(3) + else + ! Crossing into bottom mesh cell -- this is treated as incoming + ! current in (i,j,k-1) + ijk0(3) = ijk0(3) - 1 + xyz_cross(3) = xyz_cross(3) - m % width(3) + if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then + score_index = sum(t % stride(1:3) * ijk0) + IN_TOP + end if end if end if - elseif (distance == d(2)) then - if (uvw(2) > 0) then - ! Crossing into front mesh cell -- this is treated as outgoing - ! current in (i,j,k) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT - end if - ijk0(2) = ijk0(2) + 1 - xyz_cross(2) = xyz_cross(2) + m % width(2) - else - ! Crossing into back mesh cell -- this is treated as incoming - ! current in (i,j-1,k) - ijk0(2) = ijk0(2) - 1 - xyz_cross(2) = xyz_cross(2) - m % width(2) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT - end if + + ! Check for errors + if (score_index < 0 .or. score_index > t % n_total_bins) then + message = "Score index outside range." + call fatal_error() end if - else if (distance == d(3)) then - if (uvw(3) > 0) then - ! Crossing into top mesh cell -- this is treated as outgoing - ! current in (i,j,k) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP - end if - ijk0(3) = ijk0(3) + 1 - xyz_cross(3) = xyz_cross(3) + m % width(3) - else - ! Crossing into bottom mesh cell -- this is treated as incoming - ! current in (i,j,k-1) - ijk0(3) = ijk0(3) - 1 - xyz_cross(3) = xyz_cross(3) - m % width(3) - if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then - score_index = sum(t % stride(1:3) * ijk0) + IN_TOP - end if + + ! Add to surface current tally + if (score_index > 0) then + call add_to_score(t % scores(score_index, 1), p % last_wgt) end if - end if - ! Check for errors - if (score_index < 0 .or. score_index > t % n_total_bins) then - message = "Score index outside range." - call fatal_error() - end if + ! Calculate new coordinates + xyz0 = xyz0 + distance * uvw + end do - ! Add to surface current tally - if (score_index > 0) then - call add_to_score(t % scores(score_index, 1), p % last_wgt) - end if - - ! Calculate new coordinates - xyz0 = xyz0 + distance * uvw end do end subroutine score_surface_current