crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Ability to create plot.out files complete.

Browse source
This commit is contained in:
Paul Romano 2011-10-16 23:10:44 -04:00
parent 5746aa113d
commit 3a16c5875f
8 changed files with 268 additions and 63 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/constants.f90
View file

@ -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

41
src/geometry.f90
View file

@ -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)

4
src/global.f90
View file

@ -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

12
src/initialize.f90
View file

@ -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

54
src/input_xml.f90
View file

@ -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

27
src/output.f90
View file

@ -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.

186
src/plot.f90
View file

@ -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

3
src/xml-fortran/templates/plot_t.xml
View file

@ -5,6 +5,7 @@
<variable name="origin_" tag="origin" type="double-array" />
<variable name="width_" tag="width" type="double-array" />
<variable name="basis_" tag="basis" type="word" length="3" />
<variable name="basis_" tag="basis" type="word" length="3" default="'xy'" />
<variable name="pixel_" tag="pixel" type="double" default="0.01" />
</template>