diff --git a/src/constants.F90 b/src/constants.F90 index fe025c7e8a..abc0866b6c 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -32,6 +32,7 @@ module constants ! User for precision in geometry real(8), parameter :: FP_PRECISION = 1e-14_8 real(8), parameter :: FP_REL_PRECISION = 1e-5_8 + real(8), parameter :: FP_COINCIDENT = 1e-12_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 10000 diff --git a/src/geometry.F90 b/src/geometry.F90 index fecdc33cf2..fb15f6f3ce 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -47,7 +47,7 @@ contains ! Determine the specified sense of the surface in the cell and the actual ! sense of the particle with respect to the surface s => surfaces(abs(i_surface)) - actual_sense = sense(s, p % coord % xyz) + actual_sense = sense(s) specified_sense = (c % surfaces(i) > 0) ! Compare sense of point to specified sense @@ -480,14 +480,32 @@ contains ! ========================================================================== ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + ! Remove lower coordinate levels and assignment of surface + p % surface = NONE + p % coord => p % coord0 + call deallocate_coord(p % coord % next) call find_cell(found) - ! Couldn't find next cell anywhere! - if ((.not. found) .and. (run_mode /= MODE_PLOTTING)) then - message = "After particle " // trim(to_str(p % id)) // " crossed surface " & - // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & - &located in any cell and it did not leak." - call fatal_error() + 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 % coord => p % coord0 + call deallocate_coord(p % coord % next) + p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + call find_cell(found) + + ! Couldn't find next cell anywhere! This probably means there is an actual + ! undefined region in the geometry. + + if (.not. found) then + message = "After particle " // trim(to_str(p % id)) // " crossed surface " & + // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & + &located in any cell and it did not leak." + call fatal_error() + end if end if end subroutine cross_surface @@ -1141,10 +1159,9 @@ contains ! is in. !=============================================================================== - function sense(surf, xyz) result(s) + recursive function sense(surf) result(s) type(Surface), pointer :: surf ! surface - real(8), intent(in) :: xyz(3) ! coordinates of particle logical :: s ! sense of particle real(8) :: x,y,z ! coordinates of particle @@ -1163,9 +1180,9 @@ contains real(8) :: r ! radius for quadratic surfaces real(8) :: x1,y1,z1 ! upper-right corner of box - x = xyz(1) - y = xyz(2) - z = xyz(3) + x = p % coord % xyz(1) + y = p % coord % xyz(2) + z = p % coord % xyz(3) select case (surf % type) case (SURF_PX) @@ -1319,7 +1336,12 @@ contains end select ! Check which side of surface the point is on - if (func > 0) then + if (abs(func) < FP_COINCIDENT) then + ! Particle may be coincident with this surface. Artifically move the + ! particle forward a tiny bit. + p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + s = sense(surf) + elseif (func > 0) then s = .true. else s = .false. diff --git a/src/plot.F90 b/src/plot.F90 index a0ccbabe67..4e9d48c257 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -95,7 +95,7 @@ contains allocate(p) call initialize_particle() p % coord % xyz = xyz - p % coord % uvw = (/ 1, 0, 0 /) + p % coord % uvw = (/ 0.5, 0.5, 0.5 /) p % coord % universe = BASE_UNIVERSE do y = 1, img % height