diff --git a/src/constants.f90 b/src/constants.f90
index 7e26c7e9ab..225116c6d7 100644
--- a/src/constants.f90
+++ b/src/constants.f90
@@ -38,7 +38,8 @@ module constants
integer, parameter :: &
& CELL_NORMAL = 1, & ! Cell with a specified material
& CELL_FILL = 2, & ! Cell filled by a separate universe
- & CELL_LATTICE = 3 ! Cell filled with a lattice
+ & CELL_LATTICE = 3, & ! Cell filled with a lattice
+ & CELL_VOID = -1
! Lattice types
integer, parameter :: &
@@ -199,6 +200,7 @@ module constants
! Unit numbers
integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file
integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file
+ integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file
end module constants
diff --git a/src/geometry.f90 b/src/geometry.f90
index 15f6efc2b0..d48419d58b 100644
--- a/src/geometry.f90
+++ b/src/geometry.f90
@@ -200,7 +200,7 @@ contains
call message(msg)
end if
- if (surf % bc == BC_VACUUM) then
+ if (surf % bc == BC_VACUUM .and. (.not. plotting)) then
! =======================================================================
! PARTICLE LEAKS OUT OF PROBLEM
@@ -211,7 +211,7 @@ contains
end if
return
- elseif (surf % bc == BC_REFLECT) then
+ elseif (surf % bc == BC_REFLECT .and. (.not. plotting)) then
! =======================================================================
! PARTICLE REFLECTS FROM SURFACE
@@ -531,24 +531,31 @@ contains
type(Surface), pointer :: surf_p => null()
type(Lattice), pointer :: lat => null()
- cell_p => cells(p%cell)
+ if (p % cell == CELL_VOID) then
+ n_surf = n_surfaces
+ allocate(expression(n_surfaces))
+ expression = (/ (i, i=1, n_surfaces) /)
+ current_surf = 0
+ else
+ cell_p => cells(p%cell)
- current_surf = p%surface
+ current_surf = p%surface
- ! determine number of surfaces to check
- n1 = cell_p % n_surfaces
- n2 = 0
- if (cell_p % parent > 0) then
- parent_p => cells(cell_p % parent)
- n2 = parent_p % n_surfaces
- end if
- n_surf = n1 + n2
+ ! determine number of surfaces to check
+ n1 = cell_p % n_surfaces
+ n2 = 0
+ if (cell_p % parent > 0) then
+ parent_p => cells(cell_p % parent)
+ n2 = parent_p % n_surfaces
+ end if
+ n_surf = n1 + n2
- ! allocate space and assign expression
- allocate(expression(n_surf))
- expression(1:n1) = cell_p % surfaces
- if (cell_p % parent > 0) then
- expression(n1+1:n1+n2) = parent_p % surfaces
+ ! allocate space and assign expression
+ allocate(expression(n_surf))
+ expression(1:n1) = cell_p % surfaces
+ if (cell_p % parent > 0) then
+ expression(n1+1:n1+n2) = parent_p % surfaces
+ end if
end if
u = p % uvw(1)
diff --git a/src/global.f90 b/src/global.f90
index 9889ea9513..c67616f2ed 100644
--- a/src/global.f90
+++ b/src/global.f90
@@ -126,6 +126,10 @@ module global
! Plotting options
logical :: plotting = .false.
+ real(8) :: plot_origin(3)
+ real(8) :: plot_width(2)
+ real(8) :: plot_basis(6)
+ real(8) :: pixel
contains
diff --git a/src/initialize.f90 b/src/initialize.f90
index cfe4ca5151..5dacd9e56d 100644
--- a/src/initialize.f90
+++ b/src/initialize.f90
@@ -15,7 +15,7 @@ module initialize
use mcnp_random, only: RN_init_problem
use mpi_routines, only: setup_mpi
use output, only: title, echo_input, message, print_summary, &
- print_particle, header
+ print_particle, header, print_plot
use source, only: initialize_source
use string, only: int_to_str, starts_with, ends_with
use tally, only: create_tally_map, TallyObject
@@ -102,9 +102,13 @@ contains
end if
! stop timer for initialization
- if (master .and. (.not. plotting)) then
- call echo_input()
- call print_summary()
+ if (master) then
+ if (plotting) then
+ call print_plot()
+ else
+ call echo_input()
+ call print_summary()
+ end if
end if
! Stop initialization timer
diff --git a/src/input_xml.f90 b/src/input_xml.f90
index 0879edd9f5..5fa4d0b935 100644
--- a/src/input_xml.f90
+++ b/src/input_xml.f90
@@ -30,6 +30,7 @@ contains
call read_geometry_xml()
call read_materials_xml()
call read_tallies_xml()
+ if (plotting) call read_plot_xml()
end subroutine read_input_xml
@@ -782,4 +783,57 @@ contains
end subroutine read_tallies_xml
+!===============================================================================
+! READ_TALLIES_XML reads data from a tallies.xml file and parses it, checking
+! for errors and placing properly-formatted data in the right data structures
+!===============================================================================
+
+ subroutine read_plot_xml
+
+ use xml_data_plot_t
+
+ logical :: file_exists ! does tallies.xml file exist?
+ character(MAX_LINE_LEN) :: filename
+ character(MAX_LINE_LEN) :: msg
+
+ ! Check if plot.xml exists
+ filename = trim(path_input) // "plot.xml"
+ inquire(FILE=filename, EXIST=file_exists)
+ if (.not. file_exists) then
+ msg = "Plot XML file '" // trim(filename) // "' does not exist!"
+ call fatal_error(msg)
+ end if
+
+ ! Display output message
+ msg = "Reading plot XML file..."
+ call message(msg, 5)
+
+ ! Parse plot.xml file
+ call read_xml_file_plot_t(filename)
+
+ ! Copy plotting origin
+ if (size(origin_) == 3) then
+ plot_origin = origin_
+ end if
+
+ ! Copy plotting width
+ if (size(width_) == 2) then
+ plot_width = width_
+ end if
+
+ ! Read basis
+ select case (basis_)
+ case ("xy")
+ plot_basis = (/ 1, 0, 0, 0, 1, 0 /)
+ case ("yz")
+ plot_basis = (/ 0, 1, 0, 0, 0, 1 /)
+ case ("xz")
+ plot_basis = (/ 1, 0, 0, 0, 0, 1 /)
+ end select
+
+ ! Read pixel width
+ pixel = pixel_
+
+ end subroutine read_plot_xml
+
end module input_xml
diff --git a/src/output.f90 b/src/output.f90
index 800f990f6e..5fd9b5c60d 100644
--- a/src/output.f90
+++ b/src/output.f90
@@ -715,6 +715,33 @@ contains
end subroutine print_summary
+!===============================================================================
+! PRINT_PLOT displays selected options for plotting
+!===============================================================================
+
+ subroutine print_plot()
+
+ ! Display header for plotting
+ call header("PLOTTING SUMMARY")
+
+ ! Print plotting origin
+ write(ou,100) "Plotting Origin:", trim(real_to_str(plot_origin(1))) // &
+ " " // trim(real_to_str(plot_origin(2))) // " " // &
+ trim(real_to_str(plot_origin(3)))
+
+ ! Print plotting width
+ write(ou,100) "Plotting Width:", trim(real_to_str(plot_width(1))) // &
+ " " // trim(real_to_str(plot_width(2)))
+
+ ! Print pixel width
+ write(ou,100) "Pixel Width:", trim(real_to_str(pixel))
+ write(ou,*)
+
+ ! Format descriptor for columns
+100 format (1X,A,T25,A)
+
+ end subroutine print_plot
+
!===============================================================================
! PRINT_RUNTIME displays the total time elapsed for the entire run, for
! initialization, for computation, and for intercycle synchronization.
diff --git a/src/plot.f90 b/src/plot.f90
index 4add38fc42..953e7cc326 100644
--- a/src/plot.f90
+++ b/src/plot.f90
@@ -3,10 +3,10 @@ module plot
use constants
use error, only: fatal_error
use geometry, only: find_cell, dist_to_boundary, cross_surface, &
- cross_lattice
+ cross_lattice, cell_contains
use geometry_header, only: Universe, BASE_UNIVERSE
use global
- use particle_header, only: Particle
+ use particle_header, only: Particle, initialize_particle
implicit none
@@ -18,58 +18,164 @@ contains
subroutine run_plot()
- end subroutine run_plot
-
-!===============================================================================
-! TRANSPORT encompasses the main logic for moving a particle through geometry.
-!===============================================================================
-
- subroutine transport_no_collision(p)
-
- type(Particle), pointer :: p
-
- integer :: surf ! surface which particle is on
- integer :: last_cell ! most recent cell particle was in
- real(8) :: distance ! distance particle travels
- logical :: found_cell ! found cell which particle is in?
- logical :: in_lattice ! is surface crossing in lattice?
+ integer :: i ! loop index
+ integer :: surf ! surface which particle is on
+ integer :: last_cell ! most recent cell particle was in
+ real(8) :: coord(3) ! starting coordinates
+ real(8) :: last_x_coord ! bounding x coordinate
+ real(8) :: last_y_coord ! bounding y coordinate
+ real(8) :: d ! distance to boundary
+ real(8) :: distance ! distance particle travels
+ logical :: found_cell ! found cell which particle is in?
+ logical :: in_lattice ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: msg ! output/error message
- type(Universe), pointer :: univ
+ character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
+ type(Cell), pointer :: c => null()
+ type(Universe), pointer :: univ => null()
+ type(Particle), pointer :: p => null()
- if (p % cell == 0) then
+ ! Open plot file for binary writing
+ path_plot = trim(path_input) // "plot.out"
+ open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
+
+ ! Write origin, width, basis, and pixel width to file
+ write(UNIT=UNIT_PLOT) plot_origin
+ write(UNIT=UNIT_PLOT) plot_width
+ write(UNIT=UNIT_PLOT) plot_basis
+ write(UNIT=UNIT_PLOT) pixel
+
+ ! Determine coordinates of the upper-left corner of the plot
+ coord(1) = plot_origin(1) - plot_width(1) / 2.0
+ coord(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
+ coord(3) = plot_origin(3)
+
+ ! Determine bounding x and y coordinates for plot
+ last_x_coord = plot_origin(1) + plot_width(1) / 2.0
+ last_y_coord = plot_origin(2) - plot_width(2) / 2.0
+
+ ! allocate and initialize particle
+ allocate(p)
+
+ ! loop over horizontal rays
+ do while(coord(2) > last_y_coord)
+
+ ! initialize the particle and set starting coordinate and direction
+ call initialize_particle(p)
+ p % xyz = coord
+ p % xyz_local = coord
+ p % uvw = (/ 1, 0, 0 /)
+
+ ! write starting coordinate to file
+ write(UNIT=UNIT_PLOT) p % xyz
+
+ ! Find cell that particle is currently in
univ => universes(BASE_UNIVERSE)
call find_cell(univ, p, found_cell)
- ! if particle couldn't be located, print error
+ ! =======================================================================
+ ! MOVE PARTICLE FORWARD TO NEXT CELL
+
if (.not. found_cell) then
- write(msg, '(A,3ES11.3)') &
- "Could not locate cell for particle at: ", p % xyz
- call fatal_error(msg)
- end if
- end if
+ univ => universes(BASE_UNIVERSE)
+ do i = 1, univ % n_cells
+ p % xyz = coord
+ p % xyz_local = coord
+ p % cell = univ % cells(i)
- ! find energy index, interpolation factor
- do while (p % alive)
+ distance = INFINITY
+ call dist_to_boundary(p, d, surf, in_lattice)
+ if (d < distance) then
+ ! Move particle forward to next surface
+ p % xyz = p % xyz + d * p % uvw
- ! Find the distance to the nearest boundary
- call dist_to_boundary(p, distance, surf, in_lattice)
+ ! Check to make sure particle is actually going into this cell
+ ! by moving it slightly forward and seeing if the cell contains
+ ! that coordinate
- ! Advance particle
- p%xyz = p%xyz + distance * p%uvw
- p%xyz_local = p%xyz_local + distance * p%uvw
+ p % xyz = p % xyz + 1e-4 * p % uvw
+ p % xyz_local = p % xyz
- last_cell = p % cell
- p % cell = 0
- if (in_lattice) then
- p % surface = 0
- call cross_lattice(p)
- else
- p % surface = surf
+ c => cells(p % cell)
+ if (.not. cell_contains(c, p)) cycle
+
+ ! Reset coordinate to surface crossing
+ p % xyz = p % xyz - 1e-4 * p % uvw
+ p % xyz_local = p % xyz
+
+ ! Set new distance and retain pointer to this cell
+ distance = d
+ last_cell = p % cell
+ end if
+ end do
+
+ ! No cell was found on this horizontal ray
+ if (distance == INFINITY) then
+ p % xyz(1) = last_x_coord
+ p % cell = 0
+ write(UNIT_PLOT) p % xyz, p % cell
+
+ ! Move to next horizontal ray
+ coord(2) = coord(2) - pixel
+ cycle
+ end if
+
+ ! Write coordinate where next cell begins
+ write(UNIT=UNIT_PLOT) p % xyz, 0
+
+ ! Process surface crossing for next cell
+ p % cell = 0
+ p % surface = -surf
call cross_surface(p, last_cell)
end if
+
+ ! =======================================================================
+ ! MOVE PARTICLE ACROSS HORIZONTAL TRACK
+
+ do while (p % alive)
+
+ ! Calculate distance to next boundary
+ call dist_to_boundary(p, distance, surf, in_lattice)
+
+ ! Advance particle
+ p%xyz = p%xyz + distance * p%uvw
+ p%xyz_local = p%xyz_local + distance * p%uvw
+
+ ! 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
+ p % alive = .false.
+ p % xyz(1) = last_x_coord
+
+ ! If there is no cell beyond this boundary, mark it as cell 0
+ if (distance == INFINITY) p % cell = 0
+
+ ! Write ending coordinates to file
+ write(UNIT=UNIT_PLOT) p % xyz, p % cell
+ cycle
+ end if
+
+ ! Write boundary crossing coordinates to file
+ write(UNIT=UNIT_PLOT) p % xyz, p % cell
+
+ last_cell = p % cell
+ p % cell = 0
+ if (in_lattice) then
+ p % surface = 0
+ call cross_lattice(p)
+ else
+ p % surface = surf
+ call cross_surface(p, last_cell)
+ end if
+ end do
+
+ ! Move y-coordinate to next position
+ coord(2) = coord(2) - pixel
end do
- end subroutine transport_no_collision
+ ! Close plot file
+ close(UNIT=UNIT_PLOT)
+
+ end subroutine run_plot
end module plot
diff --git a/src/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plot_t.xml
index d13615fa32..b4dc9eaea8 100644
--- a/src/xml-fortran/templates/plot_t.xml
+++ b/src/xml-fortran/templates/plot_t.xml
@@ -5,6 +5,7 @@
-
+
+