mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 20:45:35 -04:00
1434 lines
45 KiB
Fortran
1434 lines
45 KiB
Fortran
module geometry
|
|
|
|
use constants
|
|
use error, only: fatal_error
|
|
use geometry_header, only: Cell, Surface, Universe, Lattice
|
|
use global
|
|
use output, only: write_message
|
|
use particle_header, only: LocalCoord, deallocate_coord
|
|
use particle_restart_write, only: write_particle_restart
|
|
use string, only: to_str
|
|
use tally, only: score_surface_current
|
|
|
|
implicit none
|
|
|
|
contains
|
|
|
|
!===============================================================================
|
|
! SIMPLE_CELL_CONTAINS determines whether a given the current coordinates of the
|
|
! particle are inside a cell defined as the intersection of a series of surfaces
|
|
!===============================================================================
|
|
|
|
function simple_cell_contains(c) result(in_cell)
|
|
|
|
type(Cell), pointer :: c
|
|
logical :: in_cell
|
|
|
|
integer :: i ! index of surfaces in cell
|
|
integer :: i_surface ! index in surfaces array (with sign)
|
|
logical :: specified_sense ! specified sense of surface in list
|
|
logical :: actual_sense ! sense of particle wrt surface
|
|
type(Surface), pointer :: s => null()
|
|
|
|
SURFACE_LOOP: do i = 1, c % n_surfaces
|
|
! Lookup surface
|
|
i_surface = c % surfaces(i)
|
|
|
|
! Check if the particle is currently on the specified surface
|
|
if (i_surface == p % surface) then
|
|
! Particle is heading into the cell
|
|
cycle
|
|
elseif (i_surface == -p % surface) then
|
|
! Particle is heading out of the cell
|
|
in_cell = .false.
|
|
return
|
|
end if
|
|
|
|
! 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)
|
|
specified_sense = (c % surfaces(i) > 0)
|
|
|
|
! Compare sense of point to specified sense
|
|
if (actual_sense .neqv. specified_sense) then
|
|
in_cell = .false.
|
|
return
|
|
end if
|
|
end do SURFACE_LOOP
|
|
|
|
! If we've reached here, then the sense matched on every surface
|
|
in_cell = .true.
|
|
|
|
end function simple_cell_contains
|
|
|
|
!===============================================================================
|
|
! FIND_CELL determines what cell a source particle is in within a particular
|
|
! universe. If the base universe is passed, the particle should be found as long
|
|
! as it's within the geometry
|
|
!===============================================================================
|
|
|
|
recursive subroutine find_cell(found, search_cells)
|
|
|
|
logical, intent(inout) :: found
|
|
integer, optional :: search_cells(:)
|
|
|
|
integer :: i ! index over cells
|
|
integer :: i_x, i_y, i_z ! indices in lattice
|
|
integer :: n_x, n_y, n_z ! size of lattice
|
|
integer :: n ! number of cells to search
|
|
integer :: index_cell ! index in cells array
|
|
real(8) :: xyz(3) ! temporary location
|
|
logical :: use_search_cells ! use cells provided as argument
|
|
type(Cell), pointer :: c ! pointer to cell
|
|
type(Lattice), pointer :: lat ! pointer to lattice
|
|
type(Universe), pointer :: univ ! universe to search in
|
|
|
|
! Remove coordinates for any lower levels
|
|
call deallocate_coord(p % coord % next)
|
|
|
|
! set size of list to search
|
|
if (present(search_cells)) then
|
|
use_search_cells = .true.
|
|
n = size(search_cells)
|
|
else
|
|
use_search_cells = .false.
|
|
univ => universes(p % coord % universe)
|
|
n = univ % n_cells
|
|
end if
|
|
|
|
do i = 1, n
|
|
! select cells based on whether we are searching a universe or a provided
|
|
! list of cells (this would be for lists of neighbor cells)
|
|
if (use_search_cells) then
|
|
index_cell = search_cells(i)
|
|
! check to make sure search cell is in same universe
|
|
if (cells(index_cell) % universe /= p % coord % universe) cycle
|
|
else
|
|
index_cell = univ % cells(i)
|
|
end if
|
|
|
|
! get pointer to cell
|
|
c => cells(index_cell)
|
|
|
|
if (simple_cell_contains(c)) then
|
|
! Set cell on this level
|
|
p % coord % cell = index_cell
|
|
|
|
! Show cell information on trace
|
|
if (verbosity >= 10 .or. trace) then
|
|
message = " Entering cell " // trim(to_str(c % id))
|
|
call write_message()
|
|
end if
|
|
|
|
if (c % type == CELL_NORMAL) then
|
|
! ====================================================================
|
|
! AT LOWEST UNIVERSE, TERMINATE SEARCH
|
|
|
|
! set material
|
|
p % last_material = p % material
|
|
p % material = c % material
|
|
|
|
elseif (c % type == CELL_FILL) then
|
|
! ====================================================================
|
|
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
|
|
|
|
! Create new level of coordinates
|
|
allocate(p % coord % next)
|
|
p % coord % next % xyz = p % coord % xyz
|
|
p % coord % next % uvw = p % coord % uvw
|
|
|
|
! Move particle to next level and set universe
|
|
p % coord => p % coord % next
|
|
p % coord % universe = c % fill
|
|
|
|
! Apply translation
|
|
if (allocated(c % translation)) then
|
|
p % coord % xyz = p % coord % xyz - c % translation
|
|
end if
|
|
|
|
! Apply rotation
|
|
if (allocated(c % rotation)) then
|
|
p % coord % xyz = matmul(c % rotation, p % coord % xyz)
|
|
p % coord % uvw = matmul(c % rotation, p % coord % uvw)
|
|
p % coord % rotated = .true.
|
|
end if
|
|
|
|
call find_cell(found)
|
|
if (.not. found) exit
|
|
|
|
elseif (c % type == CELL_LATTICE) then
|
|
! ====================================================================
|
|
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
|
|
|
|
! Set current lattice
|
|
lat => lattices(c % fill)
|
|
|
|
! determine lattice index based on position
|
|
xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
|
|
i_x = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1))
|
|
i_y = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2))
|
|
n_x = lat % dimension(1)
|
|
n_y = lat % dimension(2)
|
|
if (lat % n_dimension == 3) then
|
|
i_z = ceiling((xyz(3) - lat % lower_left(3))/lat % width(3))
|
|
n_z = lat % dimension(3)
|
|
else
|
|
i_z = 1
|
|
n_z = 1
|
|
end if
|
|
|
|
! Check if lattice coordinates are within bounds
|
|
if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
|
|
i_z < 1 .or. i_z > n_z) then
|
|
|
|
! This condition should only get hit in rare circumstances where a
|
|
! neutron hits the corner of a lattice. In this case, the neutron
|
|
! may need to be moved diagonally across the lattice. To do so, we
|
|
! remove all lower coordinate levels and then search from universe
|
|
! 0.
|
|
|
|
p % coord => p % coord0
|
|
call deallocate_coord(p % coord % next)
|
|
|
|
! Reset surface and advance particle a tiny bit
|
|
p % surface = NONE
|
|
p % coord % xyz = xyz
|
|
|
|
else
|
|
|
|
! Create new level of coordinates
|
|
allocate(p % coord % next)
|
|
|
|
! adjust local position of particle
|
|
p % coord % next % xyz(1) = p % coord % xyz(1) - &
|
|
(lat % lower_left(1) + (i_x - 0.5_8)*lat % width(1))
|
|
p % coord % next % xyz(2) = p % coord % xyz(2) - &
|
|
(lat % lower_left(2) + (i_y - 0.5_8)*lat % width(2))
|
|
if (lat % n_dimension == 3) then
|
|
p % coord % next % xyz(3) = p % coord % xyz(3) - &
|
|
(lat % lower_left(3) + (i_z - 0.5_8)*lat % width(3))
|
|
else
|
|
p % coord % next % xyz(3) = p % coord % xyz(3)
|
|
end if
|
|
p % coord % next % uvw = p % coord % uvw
|
|
|
|
! Move particle to next level
|
|
p % coord => p % coord % next
|
|
|
|
! set particle lattice indices
|
|
p % coord % lattice = c % fill
|
|
p % coord % lattice_x = i_x
|
|
p % coord % lattice_y = i_y
|
|
p % coord % lattice_z = i_z
|
|
p % coord % universe = lat % universes(i_x,i_y,i_z)
|
|
end if
|
|
|
|
call find_cell(found)
|
|
if (.not. found) exit
|
|
end if
|
|
|
|
! Found cell so we can return
|
|
found = .true.
|
|
return
|
|
end if
|
|
end do
|
|
|
|
found = .false.
|
|
|
|
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(last_cell)
|
|
|
|
integer, intent(in) :: last_cell ! last cell particle was in
|
|
|
|
real(8) :: x ! x-x0 for sphere
|
|
real(8) :: y ! y-y0 for sphere
|
|
real(8) :: z ! z-z0 for sphere
|
|
real(8) :: R ! radius of sphere
|
|
real(8) :: u ! x-component of direction
|
|
real(8) :: v ! y-component of direction
|
|
real(8) :: w ! z-component of direction
|
|
real(8) :: n1 ! x-component of surface normal
|
|
real(8) :: n2 ! y-component of surface normal
|
|
real(8) :: n3 ! z-component of surface normal
|
|
real(8) :: dot_prod ! dot product of direction and normal
|
|
real(8) :: norm ! "norm" of surface normal
|
|
integer :: i_surface ! index in surfaces
|
|
logical :: found ! particle found in universe?
|
|
type(Surface), pointer :: surf => null()
|
|
|
|
i_surface = abs(p % surface)
|
|
surf => surfaces(i_surface)
|
|
if (verbosity >= 10 .or. trace) then
|
|
message = " Crossing surface " // trim(to_str(surf % id))
|
|
call write_message()
|
|
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 % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
|
|
call score_surface_current()
|
|
end if
|
|
|
|
! Score to global leakage tally
|
|
if (tallies_on) global_tallies(LEAKAGE) % value = &
|
|
global_tallies(LEAKAGE) % value + p % wgt
|
|
|
|
! Display message
|
|
if (verbosity >= 10 .or. trace) then
|
|
message = " Leaked out of surface " // trim(to_str(surf % id))
|
|
call write_message()
|
|
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 (.not. associated(p % coord, p % coord0)) then
|
|
call write_particle_restart()
|
|
message = "Cannot reflect particle " // trim(to_str(p % id)) // &
|
|
" off surface in a lower universe."
|
|
call fatal_error()
|
|
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 in coincident with a mesh boundary
|
|
|
|
if (active_current_tallies % size() > 0) then
|
|
p % coord0 % xyz = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw
|
|
call score_surface_current()
|
|
p % coord0 % xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
|
|
end if
|
|
|
|
! Copy particle's direction cosines
|
|
u = p % coord0 % uvw(1)
|
|
v = p % coord0 % uvw(2)
|
|
w = p % coord0 % uvw(3)
|
|
|
|
select case (surf%type)
|
|
case (SURF_PX)
|
|
u = -u
|
|
|
|
case (SURF_PY)
|
|
v = -v
|
|
|
|
case (SURF_PZ)
|
|
w = -w
|
|
|
|
case (SURF_PLANE)
|
|
! Find surface coefficients and norm of vector normal to surface
|
|
n1 = surf % coeffs(1)
|
|
n2 = surf % coeffs(2)
|
|
n3 = surf % coeffs(3)
|
|
norm = n1*n1 + n2*n2 + n3*n3
|
|
dot_prod = u*n1 + v*n2 + w*n3
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*n1/norm
|
|
v = v - 2*dot_prod*n2/norm
|
|
w = w - 2*dot_prod*n3/norm
|
|
|
|
case (SURF_CYL_X)
|
|
! Find y-y0, z-z0 and dot product of direction and surface normal
|
|
y = p % coord0 % xyz(2) - surf % coeffs(1)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(2)
|
|
R = surf % coeffs(3)
|
|
dot_prod = v*y + w*z
|
|
|
|
! Reflect direction according to normal
|
|
v = v - 2*dot_prod*y/(R*R)
|
|
w = w - 2*dot_prod*z/(R*R)
|
|
|
|
case (SURF_CYL_Y)
|
|
! Find x-x0, z-z0 and dot product of direction and surface normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(2)
|
|
R = surf % coeffs(3)
|
|
dot_prod = u*x + w*z
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*x/(R*R)
|
|
w = w - 2*dot_prod*z/(R*R)
|
|
|
|
case (SURF_CYL_Z)
|
|
! Find x-x0, y-y0 and dot product of direction and surface normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
y = p % coord0 % xyz(2) - surf % coeffs(2)
|
|
R = surf % coeffs(3)
|
|
dot_prod = u*x + v*y
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*x/(R*R)
|
|
v = v - 2*dot_prod*y/(R*R)
|
|
|
|
case (SURF_SPHERE)
|
|
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
|
|
! normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
y = p % coord0 % xyz(2) - surf % coeffs(2)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(3)
|
|
R = surf % coeffs(4)
|
|
dot_prod = u*x + v*y + w*z
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*x/(R*R)
|
|
v = v - 2*dot_prod*y/(R*R)
|
|
w = w - 2*dot_prod*z/(R*R)
|
|
|
|
case (SURF_CONE_X)
|
|
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
|
|
! normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
y = p % coord0 % xyz(2) - surf % coeffs(2)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(3)
|
|
R = surf % coeffs(4)
|
|
dot_prod = (v*y + w*z - R*u*x)/((R + ONE)*R*x*x)
|
|
|
|
! Reflect direction according to normal
|
|
u = u + 2*dot_prod*R*x
|
|
v = v - 2*dot_prod*y
|
|
w = w - 2*dot_prod*z
|
|
|
|
case (SURF_CONE_Y)
|
|
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
|
|
! normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
y = p % coord0 % xyz(2) - surf % coeffs(2)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(3)
|
|
R = surf % coeffs(4)
|
|
dot_prod = (u*x + w*z - R*v*y)/((R + ONE)*R*y*y)
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*x
|
|
v = v + 2*dot_prod*R*y
|
|
w = w - 2*dot_prod*z
|
|
|
|
case (SURF_CONE_Z)
|
|
! Find x-x0, y-y0, z-z0 and dot product of direction and surface
|
|
! normal
|
|
x = p % coord0 % xyz(1) - surf % coeffs(1)
|
|
y = p % coord0 % xyz(2) - surf % coeffs(2)
|
|
z = p % coord0 % xyz(3) - surf % coeffs(3)
|
|
R = surf % coeffs(4)
|
|
dot_prod = (u*x + v*y - R*w*z)/((R + ONE)*R*z*z)
|
|
|
|
! Reflect direction according to normal
|
|
u = u - 2*dot_prod*x
|
|
v = v - 2*dot_prod*y
|
|
w = w + 2*dot_prod*R*z
|
|
|
|
case default
|
|
call write_particle_restart()
|
|
message = "Reflection not supported for surface " // &
|
|
trim(to_str(surf % id))
|
|
call fatal_error()
|
|
end select
|
|
|
|
! Set new particle direction
|
|
p % coord0 % uvw = (/ u, v, w /)
|
|
|
|
! Reassign particle's cell and surface
|
|
p % coord0 % cell = last_cell
|
|
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.
|
|
|
|
if (associated(p % coord0 % next)) then
|
|
call deallocate_coord(p % coord0 % next)
|
|
call find_cell(found)
|
|
if (.not. found) then
|
|
call write_particle_restart()
|
|
message = "Couldn't find particle after reflecting from surface."
|
|
call fatal_error()
|
|
end if
|
|
end if
|
|
|
|
! Set previous coordinate going slightly past surface crossing
|
|
p % last_xyz = p % coord0 % xyz + TINY_BIT * p % coord0 % uvw
|
|
|
|
! Diagnostic message
|
|
if (verbosity >= 10 .or. trace) then
|
|
message = " Reflected from surface " // trim(to_str(surf%id))
|
|
call write_message()
|
|
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(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(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 % coord => p % coord0
|
|
call deallocate_coord(p % coord % next)
|
|
call find_cell(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 % 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
|
|
call write_particle_restart()
|
|
message = "After particle " // trim(to_str(p % id)) // " crossed surface " &
|
|
// trim(to_str(surfaces(i_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
|
|
|
|
!===============================================================================
|
|
! CROSS_LATTICE moves a particle into a new lattice element
|
|
!===============================================================================
|
|
|
|
subroutine cross_lattice(lattice_crossed)
|
|
|
|
integer, intent(in) :: lattice_crossed
|
|
|
|
integer :: i_x, i_y, i_z ! indices in lattice
|
|
integer :: n_x, n_y, n_z ! size of lattice
|
|
real(8) :: x0, y0, z0 ! half width of lattice element
|
|
logical :: found ! particle found in cell?
|
|
type(Lattice), pointer :: lat => null()
|
|
|
|
lat => lattices(p % coord % lattice)
|
|
|
|
if (verbosity >= 10 .or. trace) then
|
|
message = " Crossing lattice " // trim(to_str(lat % id)) // &
|
|
". Current position (" // trim(to_str(p % coord % lattice_x)) &
|
|
// "," // trim(to_str(p % coord % lattice_y)) // "," // &
|
|
trim(to_str(p % coord % lattice_z)) // ")"
|
|
call write_message()
|
|
end if
|
|
|
|
if (lat % type == LATTICE_RECT) then
|
|
x0 = lat % width(1) * 0.5_8
|
|
y0 = lat % width(2) * 0.5_8
|
|
if (lat % n_dimension == 3) z0 = lat % width(3) * 0.5_8
|
|
|
|
select case (lattice_crossed)
|
|
case (LATTICE_LEFT)
|
|
! Move particle to left element
|
|
p % coord % lattice_x = p % coord % lattice_x - 1
|
|
p % coord % xyz(1) = x0
|
|
|
|
case (LATTICE_RIGHT)
|
|
! Move particle to right element
|
|
p % coord % lattice_x = p % coord % lattice_x + 1
|
|
p % coord % xyz(1) = -x0
|
|
|
|
case (LATTICE_BACK)
|
|
! Move particle to bottom element
|
|
p % coord % lattice_y = p % coord % lattice_y - 1
|
|
p % coord % xyz(2) = y0
|
|
|
|
case (LATTICE_FRONT)
|
|
! Move particle to top element
|
|
p % coord % lattice_y = p % coord % lattice_y + 1
|
|
p % coord % xyz(2) = -y0
|
|
|
|
case (LATTICE_BOTTOM)
|
|
! Move particle to bottom element
|
|
p % coord % lattice_z = p % coord % lattice_z - 1
|
|
p % coord % xyz(3) = z0
|
|
|
|
case (LATTICE_TOP)
|
|
! Move particle to top element
|
|
p % coord % lattice_z = p % coord % lattice_z + 1
|
|
p % coord % xyz(3) = -z0
|
|
|
|
end select
|
|
elseif (lat % type == LATTICE_HEX) then
|
|
! TODO: Add hex lattice support
|
|
end if
|
|
|
|
! Check to make sure still in lattice
|
|
i_x = p % coord % lattice_x
|
|
i_y = p % coord % lattice_y
|
|
i_z = p % coord % lattice_z
|
|
n_x = lat % dimension(1)
|
|
n_y = lat % dimension(2)
|
|
if (lat % n_dimension == 3) then
|
|
n_z = lat % dimension(3)
|
|
else
|
|
n_z = 1
|
|
end if
|
|
if (i_x < 1 .or. i_x > n_x .or. i_y < 1 .or. i_y > n_y .or. &
|
|
i_z < 1 .or. i_z > n_z) then
|
|
call deallocate_coord(p % coord0 % next)
|
|
p % coord => p % coord0
|
|
|
|
! Search for particle
|
|
call find_cell(found)
|
|
if (.not. found) then
|
|
call write_particle_restart()
|
|
message = "Could not locate particle " // trim(to_str(p % id)) // &
|
|
" after crossing a lattice boundary."
|
|
call fatal_error()
|
|
end if
|
|
else
|
|
! Find universe for next lattice element
|
|
p % coord % universe = lat % universes(i_x, i_y, i_z)
|
|
|
|
! Find cell in next lattice element
|
|
call find_cell(found)
|
|
if (.not. found) then
|
|
! In some circumstances, a particle crossing the corner of a cell may not
|
|
! be able to be found in the next universe. In this scenario we cut off
|
|
! all lower-level coordinates and search from universe zero
|
|
|
|
! Remove lower coordinates
|
|
call deallocate_coord(p % coord0 % next)
|
|
p % coord => p % coord0
|
|
|
|
! Search for particle
|
|
call find_cell(found)
|
|
if (.not. found) then
|
|
call write_particle_restart()
|
|
message = "Could not locate particle " // trim(to_str(p % id)) // &
|
|
" after crossing a lattice boundary."
|
|
call fatal_error()
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
end subroutine cross_lattice
|
|
|
|
!===============================================================================
|
|
! DISTANCE_TO_BOUNDARY calculates the distance to the nearest boundary for a
|
|
! particle 'p' traveling in a certain direction. For a cell in a subuniverse
|
|
! that has a parent cell, also include the surfaces of the edge of the universe.
|
|
!===============================================================================
|
|
|
|
subroutine distance_to_boundary(dist, surface_crossed, lattice_crossed)
|
|
|
|
real(8), intent(out) :: dist
|
|
integer, intent(out) :: surface_crossed
|
|
integer, intent(out) :: lattice_crossed
|
|
|
|
integer :: i ! index for surface in cell
|
|
integer :: index_surf ! index in surfaces array (with sign)
|
|
real(8) :: x,y,z ! particle coordinates
|
|
real(8) :: u,v,w ! particle directions
|
|
real(8) :: d ! evaluated distance
|
|
real(8) :: x0,y0,z0 ! coefficients for surface
|
|
real(8) :: r ! radius for quadratic surfaces
|
|
real(8) :: tmp ! dot product of surface normal with direction
|
|
real(8) :: a,b,c,k ! quadratic equation coefficients
|
|
real(8) :: quad ! discriminant of quadratic equation
|
|
logical :: on_surface ! is particle on surface?
|
|
type(Cell), pointer :: cl => null()
|
|
type(Surface), pointer :: surf => null()
|
|
type(Lattice), pointer :: lat => null()
|
|
type(LocalCoord), pointer :: coord => null()
|
|
type(LocalCoord), pointer :: final_coord => null()
|
|
|
|
! inialize distance to infinity (huge)
|
|
dist = INFINITY
|
|
lattice_crossed = NONE
|
|
nullify(final_coord)
|
|
|
|
! Get pointer to top-level coordinates
|
|
coord => p % coord0
|
|
|
|
! Loop over each universe level
|
|
LEVEL_LOOP: do while(associated(coord))
|
|
|
|
! get pointer to cell on this level
|
|
cl => cells(coord % cell)
|
|
|
|
! copy directional cosines
|
|
u = coord % uvw(1)
|
|
v = coord % uvw(2)
|
|
w = coord % uvw(3)
|
|
|
|
! =======================================================================
|
|
! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
|
|
|
|
SURFACE_LOOP: do i = 1, cl % n_surfaces
|
|
|
|
! copy local coordinates of particle
|
|
x = coord % xyz(1)
|
|
y = coord % xyz(2)
|
|
z = coord % xyz(3)
|
|
|
|
! check for coincident surface -- note that we can't skip the
|
|
! calculation in general because a particle could be on one side of a
|
|
! cylinder and still have a positive distance to the other
|
|
|
|
index_surf = cl % surfaces(i)
|
|
if (index_surf == p % surface) then
|
|
on_surface = .true.
|
|
else
|
|
on_surface = .false.
|
|
end if
|
|
|
|
! check for operators
|
|
index_surf = abs(index_surf)
|
|
if (index_surf >= OP_DIFFERENCE) cycle
|
|
|
|
! get pointer to surface
|
|
surf => surfaces(index_surf)
|
|
|
|
! TODO: Can probably combines a lot of the cases to reduce repetition
|
|
! since the algorithm is the same for (x-plane, y-plane, z-plane),
|
|
! (x-cylinder, y-cylinder, z-cylinder), etc.
|
|
|
|
select case (surf % type)
|
|
case (SURF_PX)
|
|
if (on_surface .or. u == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
x0 = surf % coeffs(1)
|
|
d = (x0 - x)/u
|
|
if (d < ZERO) d = INFINITY
|
|
end if
|
|
|
|
case (SURF_PY)
|
|
if (on_surface .or. v == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
y0 = surf % coeffs(1)
|
|
d = (y0 - y)/v
|
|
if (d < ZERO) d = INFINITY
|
|
end if
|
|
|
|
case (SURF_PZ)
|
|
if (on_surface .or. w == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
z0 = surf % coeffs(1)
|
|
d = (z0 - z)/w
|
|
if (d < ZERO) d = INFINITY
|
|
end if
|
|
|
|
case (SURF_PLANE)
|
|
A = surf % coeffs(1)
|
|
B = surf % coeffs(2)
|
|
C = surf % coeffs(3)
|
|
D = surf % coeffs(4)
|
|
|
|
tmp = A*u + B*v + C*w
|
|
if (on_surface .or. tmp == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = -(A*x + B*y + C*w - D)/tmp
|
|
if (d < ZERO) d = INFINITY
|
|
end if
|
|
|
|
case (SURF_CYL_X)
|
|
a = ONE - u*u ! v^2 + w^2
|
|
if (a == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
y0 = surf % coeffs(1)
|
|
z0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
|
|
y = y - y0
|
|
z = z - z0
|
|
k = y*v + z*w
|
|
c = y*y + z*z - r*r
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cylinder
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cylinder, thus one distance is
|
|
! positive/negative and the other is zero. The sign of k
|
|
! determines if we are facing in or out
|
|
|
|
if (k >= ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
elseif (c < ZERO) then
|
|
! particle is inside the cylinder, thus one distance must be
|
|
! negative and one must be positive. The positive distance
|
|
! will be the one with negative sign on sqrt(quad)
|
|
|
|
d = (-k + sqrt(quad))/a
|
|
|
|
else
|
|
! particle is outside the cylinder, thus both distances are
|
|
! either positive or negative. If positive, the smaller
|
|
! distance is the one with positive sign on sqrt(quad)
|
|
|
|
d = (-k - sqrt(quad))/a
|
|
if (d < ZERO) d = INFINITY
|
|
|
|
end if
|
|
end if
|
|
|
|
case (SURF_CYL_Y)
|
|
a = ONE - v*v ! u^2 + w^2
|
|
if (a == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
x0 = surf % coeffs(1)
|
|
z0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
|
|
x = x - x0
|
|
z = z - z0
|
|
k = x*u + z*w
|
|
c = x*x + z*z - r*r
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cylinder
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cylinder, thus one distance is
|
|
! positive/negative and the other is zero. The sign of k
|
|
! determines if we are facing in or out
|
|
|
|
if (k >= ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
elseif (c < ZERO) then
|
|
! particle is inside the cylinder, thus one distance must be
|
|
! negative and one must be positive. The positive distance
|
|
! will be the one with negative sign on sqrt(quad)
|
|
|
|
d = (-k + sqrt(quad))/a
|
|
|
|
else
|
|
! particle is outside the cylinder, thus both distances are
|
|
! either positive or negative. If positive, the smaller
|
|
! distance is the one with positive sign on sqrt(quad)
|
|
|
|
d = (-k - sqrt(quad))/a
|
|
if (d < ZERO) d = INFINITY
|
|
|
|
end if
|
|
end if
|
|
|
|
case (SURF_CYL_Z)
|
|
a = ONE - w*w ! u^2 + v^2
|
|
if (a == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
|
|
x = x - x0
|
|
y = y - y0
|
|
k = x*u + y*v
|
|
c = x*x + y*y - r*r
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cylinder
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cylinder, thus one distance is
|
|
! positive/negative and the other is zero. The sign of k
|
|
! determines if we are facing in or out
|
|
|
|
if (k >= ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
elseif (c < ZERO) then
|
|
! particle is inside the cylinder, thus one distance must be
|
|
! negative and one must be positive. The positive distance
|
|
! will be the one with negative sign on sqrt(quad)
|
|
|
|
d = (-k + sqrt(quad))/a
|
|
|
|
else
|
|
! particle is outside the cylinder, thus both distances are
|
|
! either positive or negative. If positive, the smaller
|
|
! distance is the one with positive sign on sqrt(quad)
|
|
|
|
d = (-k - sqrt(quad))/a
|
|
if (d <= ZERO) d = INFINITY
|
|
|
|
end if
|
|
end if
|
|
|
|
case (SURF_SPHERE)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
k = x*u + y*v + z*w
|
|
c = x*x + y*y + z*z - r*r
|
|
quad = k*k - c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with sphere
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the sphere, thus one distance is
|
|
! positive/negative and the other is zero. The sign of k
|
|
! determines if we are facing in or out
|
|
|
|
if (k >= ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = -k + sqrt(quad)
|
|
end if
|
|
|
|
elseif (c < ZERO) then
|
|
! particle is inside the sphere, thus one distance must be
|
|
! negative and one must be positive. The positive distance will
|
|
! be the one with negative sign on sqrt(quad)
|
|
|
|
d = -k + sqrt(quad)
|
|
|
|
else
|
|
! particle is outside the sphere, thus both distances are either
|
|
! positive or negative. If positive, the smaller distance is the
|
|
! one with positive sign on sqrt(quad)
|
|
|
|
d = -k - sqrt(quad)
|
|
if (d < ZERO) d = INFINITY
|
|
|
|
end if
|
|
|
|
case (SURF_CONE_X)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
a = v*v + w*w - r*u*u
|
|
k = y*v + z*w - r*x*u
|
|
c = y*y + z*z - r*x*x
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cone
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cone, thus one distance is positive/negative
|
|
! and the other is zero. The sign of k determines which distance is
|
|
! zero and which is not.
|
|
|
|
if (k >= ZERO) then
|
|
d = (-k - sqrt(quad))/a
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
else
|
|
! calculate both solutions to the quadratic
|
|
quad = sqrt(quad)
|
|
d = (-k - quad)/a
|
|
b = (-k + quad)/a
|
|
|
|
! determine the smallest positive solution
|
|
if (d < ZERO) then
|
|
if (b > ZERO) then
|
|
d = b
|
|
end if
|
|
else
|
|
if (b > ZERO) d = min(d, b)
|
|
end if
|
|
end if
|
|
|
|
! If the distance was negative, set boundary distance to infinity
|
|
if (d <= ZERO) d = INFINITY
|
|
|
|
case (SURF_CONE_Y)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
a = u*u + w*w - r*v*v
|
|
k = x*u + z*w - r*y*v
|
|
c = x*x + z*z - r*y*y
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cone
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cone, thus one distance is positive/negative
|
|
! and the other is zero. The sign of k determines which distance is
|
|
! zero and which is not.
|
|
|
|
if (k >= ZERO) then
|
|
d = (-k - sqrt(quad))/a
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
else
|
|
! calculate both solutions to the quadratic
|
|
quad = sqrt(quad)
|
|
d = (-k - quad)/a
|
|
b = (-k + quad)/a
|
|
|
|
! determine the smallest positive solution
|
|
if (d < ZERO) then
|
|
if (b > ZERO) then
|
|
d = b
|
|
end if
|
|
else
|
|
if (b > ZERO) d = min(d, b)
|
|
end if
|
|
end if
|
|
|
|
! If the distance was negative, set boundary distance to infinity
|
|
if (d <= ZERO) d = INFINITY
|
|
|
|
case (SURF_CONE_Z)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
a = u*u + v*v - r*w*w
|
|
k = x*u + y*v - r*z*w
|
|
c = x*x + y*y - r*z*z
|
|
quad = k*k - a*c
|
|
|
|
if (quad < ZERO) then
|
|
! no intersection with cone
|
|
|
|
d = INFINITY
|
|
|
|
elseif (on_surface) then
|
|
! particle is on the cone, thus one distance is positive/negative
|
|
! and the other is zero. The sign of k determines which distance is
|
|
! zero and which is not.
|
|
|
|
if (k >= ZERO) then
|
|
d = (-k - sqrt(quad))/a
|
|
else
|
|
d = (-k + sqrt(quad))/a
|
|
end if
|
|
|
|
else
|
|
! calculate both solutions to the quadratic
|
|
quad = sqrt(quad)
|
|
d = (-k - quad)/a
|
|
b = (-k + quad)/a
|
|
|
|
! determine the smallest positive solution
|
|
if (d < ZERO) then
|
|
if (b > ZERO) then
|
|
d = b
|
|
end if
|
|
else
|
|
if (b > ZERO) d = min(d, b)
|
|
end if
|
|
end if
|
|
|
|
! If the distance was negative, set boundary distance to infinity
|
|
if (d <= ZERO) d = INFINITY
|
|
|
|
end select
|
|
|
|
! Check is calculated distance is new minimum
|
|
if (d < dist) then
|
|
if (abs(d - dist)/dist >= FP_PRECISION) then
|
|
dist = d
|
|
surface_crossed = -cl % surfaces(i)
|
|
lattice_crossed = NONE
|
|
final_coord => coord
|
|
end if
|
|
end if
|
|
|
|
end do SURFACE_LOOP
|
|
|
|
! =======================================================================
|
|
! FIND MINIMUM DISTANCE TO LATTICE SURFACES
|
|
|
|
if (coord % lattice /= NONE) then
|
|
lat => lattices(coord % lattice)
|
|
if (lat % type == LATTICE_RECT) then
|
|
! copy local coordinates
|
|
x = coord % xyz(1)
|
|
y = coord % xyz(2)
|
|
z = coord % xyz(3)
|
|
|
|
! determine oncoming edge
|
|
x0 = sign(lat % width(1) * 0.5_8, u)
|
|
y0 = sign(lat % width(2) * 0.5_8, v)
|
|
|
|
! left and right sides
|
|
if (abs(x - x0) < FP_PRECISION) then
|
|
d = INFINITY
|
|
elseif (u == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (x0 - x)/u
|
|
end if
|
|
|
|
! If the lattice boundary is coincident with the parent cell boundary,
|
|
! we need to make sure that the lattice is not selected. This is
|
|
! complicated by the fact that floating point may determine that one
|
|
! is closer than the other (can't check direct equality). Thus, the
|
|
! logic here checks whether the relative difference is within floating
|
|
! point precision.
|
|
|
|
if (d < dist) then
|
|
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
|
|
dist = d
|
|
if (u > 0) then
|
|
lattice_crossed = LATTICE_RIGHT
|
|
else
|
|
lattice_crossed = LATTICE_LEFT
|
|
end if
|
|
final_coord => coord
|
|
end if
|
|
end if
|
|
|
|
! front and back sides
|
|
if (abs(y - y0) < FP_PRECISION) then
|
|
d = INFINITY
|
|
elseif (v == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (y0 - y)/v
|
|
end if
|
|
|
|
if (d < dist) then
|
|
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
|
|
dist = d
|
|
if (v > 0) then
|
|
lattice_crossed = LATTICE_FRONT
|
|
else
|
|
lattice_crossed = LATTICE_BACK
|
|
end if
|
|
final_coord => coord
|
|
end if
|
|
end if
|
|
|
|
if (lat % n_dimension == 3) then
|
|
z0 = sign(lat % width(3) * 0.5_8, w)
|
|
|
|
! top and bottom sides
|
|
if (abs(z - z0) < FP_PRECISION) then
|
|
d = INFINITY
|
|
elseif (w == ZERO) then
|
|
d = INFINITY
|
|
else
|
|
d = (z0 - z)/w
|
|
end if
|
|
|
|
if (d < dist) then
|
|
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
|
|
dist = d
|
|
if (w > 0) then
|
|
lattice_crossed = LATTICE_TOP
|
|
else
|
|
lattice_crossed = LATTICE_BOTTOM
|
|
end if
|
|
final_coord => coord
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
elseif (lat % type == LATTICE_HEX) then
|
|
! TODO: Add hex lattice support
|
|
end if
|
|
end if
|
|
|
|
coord => coord % next
|
|
|
|
end do LEVEL_LOOP
|
|
|
|
! Move particle to appropriate coordinate level
|
|
if (associated(final_coord)) p % coord => final_coord
|
|
|
|
end subroutine distance_to_boundary
|
|
|
|
!===============================================================================
|
|
! SENSE determines whether a point is on the 'positive' or 'negative' side of a
|
|
! surface. This routine is crucial for determining what cell a particular point
|
|
! is in.
|
|
!===============================================================================
|
|
|
|
recursive function sense(surf) result(s)
|
|
|
|
type(Surface), pointer :: surf ! surface
|
|
logical :: s ! sense of particle
|
|
|
|
real(8) :: x,y,z ! coordinates of particle
|
|
real(8) :: func ! surface function evaluated at point
|
|
real(8) :: A ! coefficient on x for plane
|
|
real(8) :: B ! coefficient on y for plane
|
|
real(8) :: C ! coefficient on z for plane
|
|
real(8) :: D ! coefficient for plane
|
|
real(8) :: x0,y0,z0 ! coefficients for quadratic surfaces / box
|
|
real(8) :: r ! radius for quadratic surfaces
|
|
|
|
x = p % coord % xyz(1)
|
|
y = p % coord % xyz(2)
|
|
z = p % coord % xyz(3)
|
|
|
|
select case (surf % type)
|
|
case (SURF_PX)
|
|
x0 = surf % coeffs(1)
|
|
func = x - x0
|
|
|
|
case (SURF_PY)
|
|
y0 = surf % coeffs(1)
|
|
func = y - y0
|
|
|
|
case (SURF_PZ)
|
|
z0 = surf % coeffs(1)
|
|
func = z - z0
|
|
|
|
case (SURF_PLANE)
|
|
A = surf % coeffs(1)
|
|
B = surf % coeffs(2)
|
|
C = surf % coeffs(3)
|
|
D = surf % coeffs(4)
|
|
func = A*x + B*y + C*z - D
|
|
|
|
case (SURF_CYL_X)
|
|
y0 = surf % coeffs(1)
|
|
z0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
y = y - y0
|
|
z = z - z0
|
|
func = y*y + z*z - r*r
|
|
|
|
case (SURF_CYL_Y)
|
|
x0 = surf % coeffs(1)
|
|
z0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
x = x - x0
|
|
z = z - z0
|
|
func = x*x + z*z - r*r
|
|
|
|
case (SURF_CYL_Z)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
r = surf % coeffs(3)
|
|
x = x - x0
|
|
y = y - y0
|
|
func = x*x + y*y - r*r
|
|
|
|
case (SURF_SPHERE)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
func = x*x + y*y + z*z - r*r
|
|
|
|
case (SURF_CONE_X)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
func = y*y + z*z - r*x*x
|
|
|
|
case (SURF_CONE_Y)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
func = x*x + z*z - r*y*y
|
|
|
|
case (SURF_CONE_Z)
|
|
x0 = surf % coeffs(1)
|
|
y0 = surf % coeffs(2)
|
|
z0 = surf % coeffs(3)
|
|
r = surf % coeffs(4)
|
|
x = x - x0
|
|
y = y - y0
|
|
z = z - z0
|
|
func = x*x + y*y - r*z*z
|
|
|
|
end select
|
|
|
|
! Check which side of surface the point is on
|
|
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.
|
|
end if
|
|
|
|
end function sense
|
|
|
|
!===============================================================================
|
|
! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to speed up
|
|
! searches when a cell boundary is crossed.
|
|
!===============================================================================
|
|
|
|
subroutine neighbor_lists()
|
|
|
|
integer :: i ! index in cells/surfaces array
|
|
integer :: j ! index of surface in cell
|
|
integer :: i_surface ! index in count arrays
|
|
integer, allocatable :: count_positive(:) ! # of cells on positive side
|
|
integer, allocatable :: count_negative(:) ! # of cells on negative side
|
|
logical :: positive ! positive side specified in surface list
|
|
type(Cell), pointer :: c
|
|
type(Surface), pointer :: surf
|
|
|
|
message = "Building neighboring cells lists for each surface..."
|
|
call write_message(4)
|
|
|
|
allocate(count_positive(n_surfaces))
|
|
allocate(count_negative(n_surfaces))
|
|
count_positive = 0
|
|
count_negative = 0
|
|
|
|
do i = 1, n_cells
|
|
c => cells(i)
|
|
|
|
! loop over each surface specification
|
|
do j = 1, c % n_surfaces
|
|
i_surface = c % surfaces(j)
|
|
positive = (i_surface > 0)
|
|
i_surface = abs(i_surface)
|
|
if (positive) then
|
|
count_positive(i_surface) = count_positive(i_surface) + 1
|
|
else
|
|
count_negative(i_surface) = count_negative(i_surface) + 1
|
|
end if
|
|
end do
|
|
end do
|
|
|
|
! allocate neighbor lists for each surface
|
|
do i = 1, n_surfaces
|
|
surf => surfaces(i)
|
|
if (count_positive(i) > 0) then
|
|
allocate(surf%neighbor_pos(count_positive(i)))
|
|
end if
|
|
if (count_negative(i) > 0) then
|
|
allocate(surf%neighbor_neg(count_negative(i)))
|
|
end if
|
|
end do
|
|
|
|
count_positive = 0
|
|
count_negative = 0
|
|
|
|
! loop over all cells
|
|
do i = 1, n_cells
|
|
c => cells(i)
|
|
|
|
! loop over each surface specification
|
|
do j = 1, c % n_surfaces
|
|
i_surface = c % surfaces(j)
|
|
positive = (i_surface > 0)
|
|
i_surface = abs(i_surface)
|
|
|
|
surf => surfaces(i_surface)
|
|
if (positive) then
|
|
count_positive(i_surface) = count_positive(i_surface) + 1
|
|
surf%neighbor_pos(count_positive(i_surface)) = i
|
|
else
|
|
count_negative(i_surface) = count_negative(i_surface) + 1
|
|
surf%neighbor_neg(count_negative(i_surface)) = i
|
|
end if
|
|
end do
|
|
end do
|
|
|
|
deallocate(count_positive)
|
|
deallocate(count_negative)
|
|
|
|
end subroutine neighbor_lists
|
|
|
|
end module geometry
|