Fixed geometry errors and plotting routine.

This commit is contained in:
Paul Romano 2011-12-02 11:51:52 -05:00
parent ba8169883f
commit a26ea65f2c
4 changed files with 44 additions and 28 deletions

View file

@ -101,11 +101,15 @@ contains
integer :: y ! y-index for 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.
@ -128,8 +132,6 @@ contains
! get pointer to cell
c => cells(index_cell)
! print *, 'searching cell ', c % id, cell_contains(c,p)
if (cell_contains(c, p)) then
! Set cell on this level
p % coord % cell = index_cell
@ -164,8 +166,9 @@ contains
lat => lattices(c % fill)
! determine universe based on lattice position
x = ceiling((p % coord % xyz(1) - lat % x0)/lat % width_x)
y = ceiling((p % coord % xyz(2) - lat % y0)/lat % width_y)
xyz = p % coord % xyz + TINY_BIT * p % coord % uvw
x = ceiling((xyz(1) - lat % x0)/lat % width_x)
y = ceiling((xyz(2) - lat % y0)/lat % width_y)
! Create new level of coordinates
p % in_lower_universe = .true.
@ -455,7 +458,10 @@ contains
dist = INFINITY
! left and right sides
if (u > 0) then
if (u == ZERO) then
d_left = INFINITY
d_right = INFINITY
elseif (u > 0) then
d_left = INFINITY
d_right = (x0 - x)/u
else
@ -464,7 +470,10 @@ contains
end if
! top and bottom sides
if (v > 0) then
if (v == ZERO) then
d_bottom = INFINITY
d_top = INFINITY
elseif (v > 0) then
d_bottom = INFINITY
d_top = (y0 - y)/v
else
@ -558,6 +567,7 @@ contains
! inialize distance to infinity (huge)
dist = INFINITY
lattice_crossed = .false.
! Get pointer to top-level coordinates
coord => p % coord0
@ -846,6 +856,7 @@ contains
if (d < dist) then
dist = d
surface_crossed = -cl % surfaces(i)
lattice_crossed = .false.
final_coord => coord
end if
@ -854,19 +865,22 @@ contains
! =======================================================================
! FIND MINIMUM DISTANCE TO LATTICE SURFACES
lattice_crossed = .false.
if (p % coord % lattice /= NONE) then
lat => lattices(p % coord % lattice)
if (coord % lattice /= NONE) then
lat => lattices(coord % lattice)
if (lat % type == LATTICE_RECT) then
x = p % coord % xyz(1)
y = p % coord % xyz(2)
z = p % coord % xyz(3)
! copy local coordinates
x = coord % xyz(1)
y = coord % xyz(2)
z = coord % xyz(3)
! determine oncoming edge
x0 = lat % width_x * 0.5_8
y0 = lat % width_y * 0.5_8
! left and right sides
if (u > 0) then
if (u == ZERO) then
d = INFINITY
elseif (u > 0) then
d = (x0 - x)/u
else
d = -(x + x0)/u
@ -888,7 +902,9 @@ contains
end if
! top and bottom sides
if (v > 0) then
if (v == ZERO) then
d = INFINITY
elseif (v > 0) then
d = (y0 - y)/v
else
d = -(y + y0)/v
@ -911,9 +927,6 @@ contains
end do LEVEL_LOOP
! Remove coordinates for any lower levels
call deallocate_coord(final_coord % next)
! Move particle to appropriate coordinate level
p % coord => final_coord

View file

@ -39,9 +39,9 @@ module particle_header
integer :: type ! Particle type (n, p, e, etc)
! Particle coordinates
logical :: in_lower_universe ! is particle in lower universe?
type(LocalCoord), pointer :: coord0 ! coordinates on universe 0
type(LocalCoord), pointer :: coord ! coordinates on lowest universe
logical :: in_lower_universe ! is particle in lower universe?
type(LocalCoord), pointer :: coord0 => null() ! coordinates on universe 0
type(LocalCoord), pointer :: coord => null() ! coordinates on lowest universe
! Other physical data
real(8) :: wgt ! particle weight

View file

@ -10,7 +10,7 @@ module physics
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use interpolation, only: interpolate_tab1
use output, only: write_message, print_particle
use output, only: write_message
use particle_header, only: Particle, LocalCoord
use random_lcg, only: prn
use search, only: binary_search
@ -104,7 +104,7 @@ contains
call collision(p)
! Save coordinates at collision for tallying purposes
p % last_xyz = p % coord % xyz
p % last_xyz = p % coord0 % xyz
! Set all uvws to base level -- right now, after a collision, only the
! base level uvws are changed

View file

@ -13,7 +13,9 @@ module plot
contains
!===============================================================================
! RUN_PLOT
! RUN_PLOT generates a binary stream file containing a list of surface/lattice
! crossings and what cell was traveled through. A Python script can then be used
! to generate a plot based on the recorded crossings and cells
!===============================================================================
subroutine run_plot()
@ -131,6 +133,8 @@ contains
! MOVE PARTICLE ACROSS HORIZONTAL TRACK
do while (p % alive)
! save particle's current cell
last_cell = p % coord % cell
! Calculate distance to next boundary
call distance_to_boundary(p, distance, surface_crossed, lattice_crossed)
@ -152,14 +156,13 @@ contains
if (distance == INFINITY) p % coord % cell = 0
! Write ending coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
cycle
end if
! Write boundary crossing coordinates to file
write(UNIT=UNIT_PLOT) p % coord0 % xyz, p % coord % cell
write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
last_cell = p % coord % cell
p % coord % cell = 0
if (lattice_crossed) then
p % surface = NONE
@ -170,7 +173,7 @@ contains
! Since boundary conditions are disabled in plotting mode, we need
! to manually add the last segment
if (surfaces(surface_crossed) % bc == BC_VACUUM) then
if (surfaces(abs(surface_crossed)) % bc == BC_VACUUM) then
p % coord0 % xyz(1) = last_x_coord
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
exit