diff --git a/src/geometry.F90 b/src/geometry.F90 index 1d99273de5..151a492858 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -10,7 +10,6 @@ module geometry use surface_header use stl_vector, only: VectorInt use string, only: to_str - use tally, only: score_surface_current use tally_header implicit none @@ -379,276 +378,6 @@ contains end subroutine find_cell -!=============================================================================== -! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of -! the geometry, is reflected, or crosses into a new lattice or cell -!=============================================================================== - - subroutine cross_surface(p) - type(Particle), intent(inout) :: p - - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: norm ! "norm" of surface normal - real(8) :: d ! distance between point and plane - real(8) :: xyz(3) ! Saved global coordinate - integer :: i_surface ! index in surfaces - logical :: rotational ! if rotational periodic BC applied - logical :: found ! particle found in universe? - class(Surface), pointer :: surf - - i_surface = abs(p % surface) - surf => surfaces(i_surface)%obj - if (verbosity >= 10 .or. trace) then - call write_message(" Crossing surface " // trim(to_str(surf % id))) - end if - - if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE LEAKS OUT OF PROBLEM - - ! Kill particle - p % alive = .false. - - ! Score any surface current tallies -- note that the particle is moved - ! forward slightly so that if the mesh boundary is on the surface, it is - ! still processed - - if (active_current_tallies % size() > 0) then - ! TODO: Find a better solution to score surface currents than - ! physically moving the particle forward slightly - - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - end if - - ! Score to global leakage tally - global_tally_leakage = global_tally_leakage + p % wgt - - ! Display message - if (verbosity >= 10 .or. trace) then - call write_message(" Leaked out of surface " & - &// trim(to_str(surf % id))) - end if - return - - elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE REFLECTS FROM SURFACE - - ! Do not handle reflective boundary conditions on lower universes - if (p % n_coord /= 1) then - call p % mark_as_lost("Cannot reflect particle " & - // trim(to_str(p % id)) // " off surface in a lower universe.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing is coincident with a mesh boundary - - if (active_current_tallies % size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - p % coord(1) % xyz = xyz - end if - - ! Reflect particle off surface - call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) - - ! Make sure new particle direction is normalized - u = p%coord(1)%uvw(1) - v = p%coord(1)%uvw(2) - w = p%coord(1)%uvw(3) - norm = sqrt(u*u + v*v + w*w) - p%coord(1)%uvw(:) = [u, v, w] / norm - - ! Reassign particle's cell and surface - p % coord(1) % cell = p % last_cell(p % last_n_coord) - p % surface = -p % surface - - ! If a reflective surface is coincident with a lattice or universe - ! boundary, it is necessary to redetermine the particle's coordinates in - ! the lower universes. - - p % n_coord = 1 - call find_cell(p, found) - if (.not. found) then - call p % mark_as_lost("Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf % id)) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Reflected from surface " & - &// trim(to_str(surf%id))) - end if - return - elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then - ! ======================================================================= - ! PERIODIC BOUNDARY - - ! Do not handle periodic boundary conditions on lower universes - if (p % n_coord /= 1) then - call p % mark_as_lost("Cannot transfer particle " & - // trim(to_str(p % id)) // " across surface in a lower universe.& - & Boundary conditions must be applied to universe 0.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing is coincident with a mesh boundary - - if (active_current_tallies % size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - p % coord(1) % xyz = xyz - end if - - rotational = .false. - select type (surf) - type is (SurfaceXPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceXPlane) - p % coord(1) % xyz(1) = opposite % x0 - type is (SurfaceYPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = v - p % coord(1) % uvw(2) = -u - - ! Rotate position - p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 - p % coord(1) % xyz(2) = opposite % y0 - end select - - type is (SurfaceYPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceYPlane) - p % coord(1) % xyz(2) = opposite % y0 - type is (SurfaceXPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = -v - p % coord(1) % uvw(2) = u - - ! Rotate position - p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 - p % coord(1) % xyz(1) = opposite % x0 - end select - - type is (SurfaceZPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceZPlane) - p % coord(1) % xyz(3) = opposite % z0 - end select - - type is (SurfacePlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfacePlane) - ! Get surface normal for opposite plane - xyz(:) = opposite % normal(p % coord(1) % xyz) - - ! Determine distance to plane - norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3) - d = opposite % evaluate(p % coord(1) % xyz) / norm - - ! Move particle along normal vector based on distance - p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz - end select - end select - - ! Reassign particle's surface - if (rotational) then - p % surface = surf % i_periodic - else - p % surface = sign(surf % i_periodic, p % surface) - end if - - ! Figure out what cell particle is in now - p % n_coord = 1 - call find_cell(p, found) - if (.not. found) then - call p % mark_as_lost("Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf % id)) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Hit periodic boundary on surface " & - // trim(to_str(surf%id))) - end if - return - end if - - ! ========================================================================== - ! SEARCH NEIGHBOR LISTS FOR NEXT CELL - - if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then - ! If coming from negative side of surface, search all the neighboring - ! cells on the positive side - - call find_cell(p, found, surf%neighbor_pos) - if (found) return - - elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then - ! If coming from positive side of surface, search all the neighboring - ! cells on the negative side - - call find_cell(p, found, surf%neighbor_neg) - if (found) return - - end if - - ! ========================================================================== - ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS - - ! Remove lower coordinate levels and assignment of surface - p % surface = NONE - p % n_coord = 1 - call find_cell(p, found) - - if (run_mode /= MODE_PLOTTING .and. (.not. found)) then - ! If a cell is still not found, there are two possible causes: 1) there is - ! a void in the model, and 2) the particle hit a surface at a tangent. If - ! the particle is really traveling tangent to a surface, if we move it - ! forward a tiny bit it should fix the problem. - - p % n_coord = 1 - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call find_cell(p, found) - - ! Couldn't find next cell anywhere! This probably means there is an actual - ! undefined region in the geometry. - - if (.not. found) then - call p % mark_as_lost("After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf % id)) & - // " it could not be located in any cell and it did not leak.") - return - end if - end if - - end subroutine cross_surface - !=============================================================================== ! CROSS_LATTICE moves a particle into a new lattice element !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index 726a263adf..d110af4622 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -4,8 +4,8 @@ module tracking use cross_section, only: calculate_xs use error, only: fatal_error, warning use geometry_header, only: cells - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice, check_cell_overlap + use geometry, only: find_cell, distance_to_boundary, cross_lattice, & + check_cell_overlap use output, only: write_message use message_passing use mgxs_header @@ -17,6 +17,7 @@ module tracking use settings use simulation_header use string, only: to_str + use surface_header use tally_header use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & @@ -278,4 +279,274 @@ contains end subroutine transport +!=============================================================================== +! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of +! the geometry, is reflected, or crosses into a new lattice or cell +!=============================================================================== + + subroutine cross_surface(p) + type(Particle), intent(inout) :: p + + real(8) :: u ! x-component of direction + real(8) :: v ! y-component of direction + real(8) :: w ! z-component of direction + real(8) :: norm ! "norm" of surface normal + real(8) :: d ! distance between point and plane + real(8) :: xyz(3) ! Saved global coordinate + integer :: i_surface ! index in surfaces + logical :: rotational ! if rotational periodic BC applied + logical :: found ! particle found in universe? + class(Surface), pointer :: surf + + i_surface = abs(p % surface) + surf => surfaces(i_surface)%obj + if (verbosity >= 10 .or. trace) then + call write_message(" Crossing surface " // trim(to_str(surf % id))) + end if + + if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then + ! ======================================================================= + ! PARTICLE LEAKS OUT OF PROBLEM + + ! Kill particle + p % alive = .false. + + ! Score any surface current tallies -- note that the particle is moved + ! forward slightly so that if the mesh boundary is on the surface, it is + ! still processed + + if (active_current_tallies % size() > 0) then + ! TODO: Find a better solution to score surface currents than + ! physically moving the particle forward slightly + + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + end if + + ! Score to global leakage tally + global_tally_leakage = global_tally_leakage + p % wgt + + ! Display message + if (verbosity >= 10 .or. trace) then + call write_message(" Leaked out of surface " & + &// trim(to_str(surf % id))) + end if + return + + elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then + ! ======================================================================= + ! PARTICLE REFLECTS FROM SURFACE + + ! Do not handle reflective boundary conditions on lower universes + if (p % n_coord /= 1) then + call p % mark_as_lost("Cannot reflect particle " & + // trim(to_str(p % id)) // " off surface in a lower universe.") + return + end if + + ! Score surface currents since reflection causes the direction of the + ! particle to change -- artificially move the particle slightly back in + ! case the surface crossing is coincident with a mesh boundary + + if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + p % coord(1) % xyz = xyz + end if + + ! Reflect particle off surface + call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) + + ! Make sure new particle direction is normalized + u = p%coord(1)%uvw(1) + v = p%coord(1)%uvw(2) + w = p%coord(1)%uvw(3) + norm = sqrt(u*u + v*v + w*w) + p%coord(1)%uvw(:) = [u, v, w] / norm + + ! Reassign particle's cell and surface + p % coord(1) % cell = p % last_cell(p % last_n_coord) + p % surface = -p % surface + + ! If a reflective surface is coincident with a lattice or universe + ! boundary, it is necessary to redetermine the particle's coordinates in + ! the lower universes. + + p % n_coord = 1 + call find_cell(p, found) + if (.not. found) then + call p % mark_as_lost("Couldn't find particle after reflecting& + & from surface " // trim(to_str(surf % id)) // ".") + return + end if + + ! Set previous coordinate going slightly past surface crossing + p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + + ! Diagnostic message + if (verbosity >= 10 .or. trace) then + call write_message(" Reflected from surface " & + &// trim(to_str(surf%id))) + end if + return + elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then + ! ======================================================================= + ! PERIODIC BOUNDARY + + ! Do not handle periodic boundary conditions on lower universes + if (p % n_coord /= 1) then + call p % mark_as_lost("Cannot transfer particle " & + // trim(to_str(p % id)) // " across surface in a lower universe.& + & Boundary conditions must be applied to universe 0.") + return + end if + + ! Score surface currents since reflection causes the direction of the + ! particle to change -- artificially move the particle slightly back in + ! case the surface crossing is coincident with a mesh boundary + + if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + p % coord(1) % xyz = xyz + end if + + rotational = .false. + select type (surf) + type is (SurfaceXPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceXPlane) + p % coord(1) % xyz(1) = opposite % x0 + type is (SurfaceYPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = v + p % coord(1) % uvw(2) = -u + + ! Rotate position + p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 + p % coord(1) % xyz(2) = opposite % y0 + end select + + type is (SurfaceYPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceYPlane) + p % coord(1) % xyz(2) = opposite % y0 + type is (SurfaceXPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = -v + p % coord(1) % uvw(2) = u + + ! Rotate position + p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 + p % coord(1) % xyz(1) = opposite % x0 + end select + + type is (SurfaceZPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceZPlane) + p % coord(1) % xyz(3) = opposite % z0 + end select + + type is (SurfacePlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfacePlane) + ! Get surface normal for opposite plane + xyz(:) = opposite % normal(p % coord(1) % xyz) + + ! Determine distance to plane + norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3) + d = opposite % evaluate(p % coord(1) % xyz) / norm + + ! Move particle along normal vector based on distance + p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz + end select + end select + + ! Reassign particle's surface + if (rotational) then + p % surface = surf % i_periodic + else + p % surface = sign(surf % i_periodic, p % surface) + end if + + ! Figure out what cell particle is in now + p % n_coord = 1 + call find_cell(p, found) + if (.not. found) then + call p % mark_as_lost("Couldn't find particle after hitting & + &periodic boundary on surface " // trim(to_str(surf % id)) // ".") + return + end if + + ! Set previous coordinate going slightly past surface crossing + p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + + ! Diagnostic message + if (verbosity >= 10 .or. trace) then + call write_message(" Hit periodic boundary on surface " & + // trim(to_str(surf%id))) + end if + return + end if + + ! ========================================================================== + ! SEARCH NEIGHBOR LISTS FOR NEXT CELL + + if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then + ! If coming from negative side of surface, search all the neighboring + ! cells on the positive side + + call find_cell(p, found, surf%neighbor_pos) + if (found) return + + elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then + ! If coming from positive side of surface, search all the neighboring + ! cells on the negative side + + call find_cell(p, found, surf%neighbor_neg) + if (found) return + + end if + + ! ========================================================================== + ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + + ! Remove lower coordinate levels and assignment of surface + p % surface = NONE + p % n_coord = 1 + call find_cell(p, found) + + if (run_mode /= MODE_PLOTTING .and. (.not. found)) then + ! If a cell is still not found, there are two possible causes: 1) there is + ! a void in the model, and 2) the particle hit a surface at a tangent. If + ! the particle is really traveling tangent to a surface, if we move it + ! forward a tiny bit it should fix the problem. + + p % n_coord = 1 + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + call find_cell(p, found) + + ! Couldn't find next cell anywhere! This probably means there is an actual + ! undefined region in the geometry. + + if (.not. found) then + call p % mark_as_lost("After particle " // trim(to_str(p % id)) & + // " crossed surface " // trim(to_str(surf % id)) & + // " it could not be located in any cell and it did not leak.") + return + end if + end if + + end subroutine cross_surface + end module tracking