Fixed cross_surface to handle lattices when searching neighbor lists. Geometry plotter can now handle lattices.

This commit is contained in:
Paul Romano 2011-10-17 13:10:41 -04:00
parent a4d544671c
commit 417d36b9c7
2 changed files with 65 additions and 4 deletions

View file

@ -176,6 +176,8 @@ contains
integer :: i ! index of neighbors
integer :: index_cell ! index in cells array
integer :: i_x ! x index in lattice
integer :: i_y ! y index in lattice
real(8) :: x ! x-x0 for sphere
real(8) :: y ! y-y0 for sphere
real(8) :: z ! z-z0 for sphere
@ -190,8 +192,9 @@ contains
real(8) :: norm ! "norm" of surface normal
logical :: found ! particle found in universe?
character(MAX_LINE_LEN) :: msg ! output/error message?
type(Surface), pointer :: surf
type(Cell), pointer :: c
type(Surface), pointer :: surf => null()
type(Cell), pointer :: c => null()
type(Lattice), pointer :: lat => null()
type(Universe), pointer :: lower_univ => null()
surf => surfaces(abs(p % surface))
@ -332,6 +335,31 @@ contains
msg = "Could not locate particle in universe: "
call fatal_error(msg)
end if
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill)
p % lattice = c % fill
! determine universe based on lattice position
i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x)
i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y)
lower_univ => universes(lat % element(i_x,i_y))
! adjust local position of particle
p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x)
p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y)
p%xyz_local(3) = p%xyz(3)
! set particle lattice indices
p % index_x = i_x
p % index_y = i_y
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call fatal_error(msg)
end if
else
! set current pointers
p % cell = index_cell
@ -354,6 +382,31 @@ contains
msg = "Could not locate particle in universe: "
call fatal_error(msg)
end if
elseif (c % type == CELL_LATTICE) then
! Set current lattice
lat => lattices(c % fill)
p % lattice = c % fill
! determine universe based on lattice position
i_x = ceiling((p%xyz(1) - lat%x0)/lat%width_x)
i_y = ceiling((p%xyz(2) - lat%y0)/lat%width_y)
lower_univ => universes(lat % element(i_x,i_y))
! adjust local position of particle
p%xyz_local(1) = p%xyz(1) - (lat%x0 + (i_x-0.5)*lat%width_x)
p%xyz_local(2) = p%xyz(2) - (lat%y0 + (i_y-0.5)*lat%width_y)
p%xyz_local(3) = p%xyz(3)
! set particle lattice indices
p % index_x = i_x
p % index_y = i_y
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call fatal_error(msg)
end if
else
! set current pointers
p % cell = index_cell

View file

@ -142,7 +142,7 @@ contains
! If next boundary crossing is out of range of the plot, only include
! the visible portion and move to next horizontal ray
if (p % xyz(1) > last_x_coord) then
if (p % xyz(1) >= last_x_coord) then
p % alive = .false.
p % xyz(1) = last_x_coord
@ -165,8 +165,16 @@ contains
else
p % surface = surf
call cross_surface(p, last_cell)
! Since boundary conditions are disabled in plotting mode, we need
! to manually add the last segment
if (surfaces(surf) % bc == BC_VACUUM) then
p % xyz(1) = last_x_coord
write(UNIT=UNIT_PLOT) p % xyz, 0
exit
end if
end if
end do
! Move y-coordinate to next position