mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
finalized ppm plotting
This commit is contained in:
nhorelik
parent
6a72d02d8d
commit
114aebc487
13 changed files with 429 additions and 194 deletions
|
|
@ -68,6 +68,7 @@ global.o: geometry_header.o
|
|||
global.o: material_header.o
|
||||
global.o: mesh_header.o
|
||||
global.o: particle_header.o
|
||||
global.o: plot_header.o
|
||||
global.o: source_header.o
|
||||
global.o: tally_header.o
|
||||
global.o: timing.o
|
||||
|
|
@ -106,6 +107,7 @@ input_xml.o: geometry_header.o
|
|||
input_xml.o: global.o
|
||||
input_xml.o: mesh_header.o
|
||||
input_xml.o: output.o
|
||||
input_xml.o: plot_header.o
|
||||
input_xml.o: string.o
|
||||
input_xml.o: tally_header.o
|
||||
input_xml.o: xml-fortran/templates/cross_sections_t.o
|
||||
|
|
@ -155,6 +157,7 @@ output.o: geometry_header.o
|
|||
output.o: global.o
|
||||
output.o: mesh_header.o
|
||||
output.o: particle_header.o
|
||||
output.o: plot_header.o
|
||||
output.o: string.o
|
||||
output.o: tally_header.o
|
||||
|
||||
|
|
@ -182,8 +185,12 @@ plot.o: error.o
|
|||
plot.o: geometry.o
|
||||
plot.o: geometry_header.o
|
||||
plot.o: global.o
|
||||
plot.o: output.o
|
||||
plot.o: particle_header.o
|
||||
plot.o: plot_header.o
|
||||
plot.o: ppmlib.o
|
||||
plot.o: source.o
|
||||
plot.o: string.o
|
||||
|
||||
search.o: constants.o
|
||||
search.o: error.o
|
||||
|
|
|
|||
|
|
@ -28,7 +28,7 @@ USE_HDF5 = no
|
|||
|
||||
ifeq ($(COMPILER),gnu)
|
||||
F90 = gfortran
|
||||
F90FLAGS := -cpp -fbacktrace
|
||||
F90FLAGS := -cpp -fbacktrace -DNO_F2008
|
||||
LDFLAGS =
|
||||
|
||||
# Debugging
|
||||
|
|
|
|||
|
|
@ -28,6 +28,8 @@ output.o \
|
|||
particle_header.o \
|
||||
physics.o \
|
||||
plot.o \
|
||||
plot_header.o \
|
||||
ppmlib.o \
|
||||
random_lcg.o \
|
||||
search.o \
|
||||
source.o \
|
||||
|
|
@ -35,4 +37,4 @@ source_header.o \
|
|||
string.o \
|
||||
tally.o \
|
||||
tally_header.o \
|
||||
timing.o
|
||||
timing.o
|
||||
|
|
|
|||
|
|
@ -64,6 +64,7 @@ module geometry_header
|
|||
& surfaces(:) ! List of surfaces bounding cell -- note that
|
||||
! parentheses, union, etc operators will be listed
|
||||
! here too
|
||||
integer, allocatable :: rgb(:) ! plotting color
|
||||
end type Cell
|
||||
|
||||
! array index of universe 0
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@ module global
|
|||
use material_header, only: Material
|
||||
use mesh_header, only: StructuredMesh
|
||||
use particle_header, only: Particle
|
||||
use plot_header, only: Plot
|
||||
use source_header, only: ExtSource
|
||||
use tally_header, only: TallyObject, TallyMap, TallyScore
|
||||
use timing, only: Timer
|
||||
|
|
@ -38,6 +39,7 @@ module global
|
|||
type(Lattice), allocatable, target :: lattices(:)
|
||||
type(Surface), allocatable, target :: surfaces(:)
|
||||
type(Material), allocatable, target :: materials(:)
|
||||
type(Plot), allocatable, target :: plots(:)
|
||||
|
||||
! Size of main arrays
|
||||
integer :: n_cells ! # of cells
|
||||
|
|
@ -45,6 +47,7 @@ module global
|
|||
integer :: n_lattices ! # of lattices
|
||||
integer :: n_surfaces ! # of surfaces
|
||||
integer :: n_materials ! # of materials
|
||||
integer :: n_plots ! # of plots
|
||||
|
||||
! These dictionaries provide a fast lookup mechanism -- the key is the
|
||||
! user-specified identifier and the value is the index in the corresponding
|
||||
|
|
@ -187,10 +190,6 @@ module global
|
|||
! PLOTTING VARIABLES
|
||||
|
||||
logical :: plotting = .false.
|
||||
real(8) :: plot_origin(3)
|
||||
real(8) :: plot_width(2)
|
||||
real(8) :: plot_basis(6)
|
||||
real(8) :: pixel
|
||||
|
||||
! ============================================================================
|
||||
! HDF5 VARIABLES
|
||||
|
|
@ -235,6 +234,7 @@ contains
|
|||
if (allocated(lattices)) deallocate(lattices)
|
||||
if (allocated(surfaces)) deallocate(surfaces)
|
||||
if (allocated(materials)) deallocate(materials)
|
||||
if (allocated(plots)) deallocate(plots)
|
||||
|
||||
! Deallocate cross section data, listings, and cache
|
||||
if (allocated(nuclides)) deallocate(nuclides)
|
||||
|
|
|
|||
|
|
@ -8,6 +8,8 @@ module input_xml
|
|||
use global
|
||||
use mesh_header, only: StructuredMesh
|
||||
use output, only: write_message
|
||||
use plot_header
|
||||
use random_lcg, only: prn
|
||||
use string, only: lower_case, to_str, str_to_int, str_to_real, &
|
||||
split_string, starts_with, ends_with
|
||||
use tally_header, only: TallyObject
|
||||
|
|
@ -259,6 +261,14 @@ contains
|
|||
c % material = cell_(i) % material
|
||||
c % fill = cell_(i) % fill
|
||||
|
||||
! Set plot color
|
||||
if (plotting) then
|
||||
allocate(c % rgb(3))
|
||||
c % rgb(1) = prn()*255
|
||||
c % rgb(2) = prn()*255
|
||||
c % rgb(3) = prn()*255
|
||||
end if
|
||||
|
||||
! Check to make sure that either material or fill was specified
|
||||
if (c % material == 0 .and. c % fill == 0) then
|
||||
message = "Neither material nor fill was specified for cell " // &
|
||||
|
|
@ -529,6 +539,14 @@ contains
|
|||
! Copy material id
|
||||
m % id = material_(i) % id
|
||||
|
||||
! Set plot color
|
||||
if (plotting) then
|
||||
allocate(m % rgb(3))
|
||||
m % rgb(1) = prn()*255
|
||||
m % rgb(2) = prn()*255
|
||||
m % rgb(3) = prn()*255
|
||||
end if
|
||||
|
||||
! Copy density -- the default value for the units is given in the
|
||||
! material_t.xml file and doesn't need to be specified here, hence case
|
||||
! default results in an error.
|
||||
|
|
@ -1086,16 +1104,17 @@ 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
|
||||
! READ_PLOT_XML reads data from a plot.xml file
|
||||
!===============================================================================
|
||||
|
||||
subroutine read_plot_xml
|
||||
|
||||
use xml_data_plot_t
|
||||
|
||||
logical :: file_exists ! does tallies.xml file exist?
|
||||
character(MAX_LINE_LEN) :: filename ! absolute path to tallies.xml
|
||||
integer i
|
||||
logical :: file_exists ! does plot.xml file exist?
|
||||
character(MAX_LINE_LEN) :: filename ! absolute path to plot.xml
|
||||
type(Plot), pointer :: pl => null()
|
||||
|
||||
! Check if plot.xml exists
|
||||
filename = trim(path_input) // "plot.xml"
|
||||
|
|
@ -1112,28 +1131,79 @@ contains
|
|||
! Parse plot.xml file
|
||||
call read_xml_file_plot_t(filename)
|
||||
|
||||
! Copy plotting origin
|
||||
if (size(origin_) == 3) then
|
||||
plot_origin = origin_
|
||||
end if
|
||||
! Allocate plots array
|
||||
n_plots = size(plot_)
|
||||
allocate(plots(n_plots))
|
||||
|
||||
! Copy plotting width
|
||||
if (size(width_) == 2) then
|
||||
plot_width = width_
|
||||
end if
|
||||
do i = 1, n_plots
|
||||
pl => plots(i)
|
||||
|
||||
! 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
|
||||
! Copy data into plots
|
||||
pl % id = plot_(i) % id
|
||||
pl % aspect = plot_(i) % aspect
|
||||
|
||||
! Read pixel width
|
||||
pixel = pixel_
|
||||
if (size(plot_(i) % pixels) == 2) then
|
||||
pl % pixels = plot_(i) % pixels
|
||||
end if
|
||||
|
||||
select case (plot_(i) % color)
|
||||
case ("cell")
|
||||
pl % color = PLOT_COLOR_CELLS
|
||||
case ("mat")
|
||||
pl % color = PLOT_COLOR_MATS
|
||||
case default
|
||||
message = "Unsupported plot color '" // plot_(i) % color &
|
||||
// "' in plot " // trim(to_str(i))
|
||||
call fatal_error()
|
||||
end select
|
||||
|
||||
select case (plot_(i) % type)
|
||||
case ("slice")
|
||||
pl % type = PLOT_TYPE_SLICE
|
||||
!case ("points")
|
||||
! pl % type = PLOT_TYPE_POINTS
|
||||
case default
|
||||
message = "Unsupported plot type '" // plot_(i) % type &
|
||||
// "' in plot " // trim(to_str(i))
|
||||
call fatal_error()
|
||||
end select
|
||||
|
||||
select case (plot_(i) % basis)
|
||||
case ("xy")
|
||||
pl % basis = PLOT_BASIS_XY
|
||||
case ("xz")
|
||||
pl % basis = PLOT_BASIS_XZ
|
||||
case ("yz")
|
||||
pl % basis = PLOT_BASIS_YZ
|
||||
case default
|
||||
message = "Unsupported plot basis '" // plot_(i) % basis &
|
||||
// "' in plot " // trim(to_str(i))
|
||||
call fatal_error()
|
||||
end select
|
||||
|
||||
! Copy plotting origin
|
||||
if (size(plot_(i) % origin) == 3) then
|
||||
pl % origin = plot_(i) % origin
|
||||
else
|
||||
message = "Origin must be length 3 " &
|
||||
// "in plot " // trim(to_str(i))
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Copy plotting width
|
||||
if (size(plot_(i) % width) == 3) then
|
||||
pl % width = plot_(i) % width
|
||||
else if (size(plot_(i) % width) == 2) then
|
||||
pl % width(1) = plot_(i) % width(1)
|
||||
pl % width(2) = plot_(i) % width(2)
|
||||
else
|
||||
message = "Bad plot width " &
|
||||
// "in plot " // trim(to_str(i))
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
|
||||
end do
|
||||
|
||||
end subroutine read_plot_xml
|
||||
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@ program main
|
|||
use initialize, only: initialize_run
|
||||
use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
|
||||
use output, only: write_message, header
|
||||
use plot, only: run_plot
|
||||
use plotter, only: run_plot
|
||||
use physics, only: transport
|
||||
use random_lcg, only: set_particle_seed
|
||||
use source, only: get_source_particle
|
||||
|
|
|
|||
|
|
@ -21,6 +21,8 @@ module material_header
|
|||
character(12) :: sab_name ! name of S(a,b) table
|
||||
integer :: sab_table = 0 ! index in sab_tables
|
||||
integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table
|
||||
|
||||
integer, allocatable :: rgb(:) ! plotting color
|
||||
end type Material
|
||||
|
||||
end module material_header
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@ module output
|
|||
use global
|
||||
use mesh_header, only: StructuredMesh
|
||||
use particle_header, only: LocalCoord
|
||||
use plot_header
|
||||
use string, only: upper_case, to_str
|
||||
use tally_header, only: TallyObject
|
||||
|
||||
|
|
@ -891,21 +892,46 @@ contains
|
|||
|
||||
subroutine print_plot()
|
||||
|
||||
integer i
|
||||
type(Plot), pointer :: pl => null()
|
||||
|
||||
! Display header for plotting
|
||||
call header("PLOTTING SUMMARY")
|
||||
|
||||
! Print plotting origin
|
||||
write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // &
|
||||
" " // trim(to_str(plot_origin(2))) // " " // &
|
||||
trim(to_str(plot_origin(3)))
|
||||
do i=1,n_plots
|
||||
pl => plots(i)
|
||||
|
||||
! Print plotting width
|
||||
write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // &
|
||||
" " // trim(to_str(plot_width(2)))
|
||||
! Print plot id
|
||||
write(ou,100) "Plot ID:", trim(to_str(pl % id))
|
||||
|
||||
! Print pixel width
|
||||
write(ou,100) "Pixel Width:", trim(to_str(pixel))
|
||||
write(ou,*)
|
||||
! Print plotting origin
|
||||
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
|
||||
" " // trim(to_str(pl % origin(2))) // " " // &
|
||||
trim(to_str(pl % origin(3)))
|
||||
|
||||
! Print plotting width
|
||||
if (pl % type == PLOT_TYPE_SLICE) then
|
||||
|
||||
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
|
||||
" " // trim(to_str(pl % width(2)))
|
||||
write(ou,100) "Coloring:", trim(to_str(pl % color))
|
||||
write(ou,100) "Basis:", trim(to_str(pl % basis))
|
||||
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
|
||||
trim(to_str(pl % pixels(2)))
|
||||
|
||||
else if (pl % type == PLOT_TYPE_POINTS) then
|
||||
|
||||
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
|
||||
" " // trim(to_str(pl % width(2))) // " " &
|
||||
// trim(to_str(pl % width(3)))
|
||||
write(ou,100) "Coloring:", trim(to_str(pl % color))
|
||||
write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect))
|
||||
|
||||
end if
|
||||
|
||||
write(ou,*)
|
||||
|
||||
end do
|
||||
|
||||
! Format descriptor for columns
|
||||
100 format (1X,A,T25,A)
|
||||
|
|
|
|||
287
src/plot.F90
287
src/plot.F90
|
|
@ -1,4 +1,4 @@
|
|||
module plot
|
||||
module plotter
|
||||
|
||||
use constants
|
||||
use error, only: fatal_error
|
||||
|
|
@ -6,185 +6,172 @@ module plot
|
|||
cross_lattice, cell_contains
|
||||
use geometry_header, only: Universe, BASE_UNIVERSE
|
||||
use global
|
||||
use output, only: write_message
|
||||
use particle_header, only: LocalCoord, deallocate_coord
|
||||
use plot_header
|
||||
use ppmlib, only: Image, init_image, allocate_image, &
|
||||
deallocate_image, set_pixel
|
||||
use source, only: initialize_particle
|
||||
use string, only: to_str
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
|
||||
!===============================================================================
|
||||
! 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
|
||||
! RUN_PLOT
|
||||
!===============================================================================
|
||||
|
||||
subroutine run_plot()
|
||||
|
||||
integer :: i ! loop index
|
||||
integer :: surface_crossed ! surface which particle is on
|
||||
integer :: lattice_crossed ! is surface crossing in lattice?
|
||||
integer :: last_cell ! most recent cell particle was in
|
||||
integer :: enter_surface ! entrance surface
|
||||
real(8) :: xyz(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?
|
||||
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
|
||||
type(Cell), pointer :: c => null()
|
||||
type(Universe), pointer :: univ => null()
|
||||
type(LocalCoord), pointer :: coord => null()
|
||||
integer :: i
|
||||
type(Plot), pointer :: pl => null()
|
||||
|
||||
! Open plot file for binary writing
|
||||
path_plot = trim(path_input) // "plot.out"
|
||||
open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
|
||||
do i=1,n_plots
|
||||
pl => plots(i)
|
||||
|
||||
! 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
|
||||
! Display output message
|
||||
message = "Processing plot " // trim(to_str(pl % id)) // "..."
|
||||
call write_message(5)
|
||||
|
||||
! Determine coordinates of the upper-left corner of the plot
|
||||
xyz(1) = plot_origin(1) - plot_width(1) / 2.0
|
||||
xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
|
||||
xyz(3) = plot_origin(3)
|
||||
if (pl % type == PLOT_TYPE_SLICE) then
|
||||
|
||||
! 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
|
||||
! create 2d image
|
||||
call create_ppm(pl)
|
||||
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine run_plot
|
||||
|
||||
!===============================================================================
|
||||
! create_ppm
|
||||
!===============================================================================
|
||||
subroutine create_ppm(pl)
|
||||
|
||||
type(Plot), pointer :: pl
|
||||
|
||||
type(Image) :: img
|
||||
integer :: in_i, out_i
|
||||
integer :: x, y
|
||||
integer :: r, g, b
|
||||
real(8) :: in_pixel, out_pixel
|
||||
real(8) :: xyz(3)
|
||||
logical :: found_cell
|
||||
type(Cell), pointer :: c => null()
|
||||
|
||||
call init_image(img)
|
||||
call allocate_image(img, pl % pixels(1), pl % pixels(2))
|
||||
|
||||
if (pl % basis == PLOT_BASIS_XY) then
|
||||
in_i = 1
|
||||
out_i = 2
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(2)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
|
||||
xyz(2) = pl % origin(2) - pl % width(2) / 2.0
|
||||
xyz(3) = pl % origin(3)
|
||||
else if (pl % basis == PLOT_BASIS_XZ) then
|
||||
in_i = 1
|
||||
out_i = 3
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(3)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
|
||||
xyz(2) = pl % origin(2)
|
||||
xyz(3) = pl % origin(3) - pl % width(3) / 2.0
|
||||
else if (pl % basis == PLOT_BASIS_YZ) then
|
||||
in_i = 2
|
||||
out_i = 3
|
||||
in_pixel = pl % width(2)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(3)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1)
|
||||
xyz(2) = pl % origin(2) - pl % width(2) / 2.0
|
||||
xyz(3) = pl % origin(3) - pl % width(3) / 2.0
|
||||
end if
|
||||
|
||||
! allocate and initialize particle
|
||||
allocate(p)
|
||||
call initialize_particle()
|
||||
p % coord % xyz = xyz
|
||||
p % coord % uvw = (/ 1, 0, 0 /)
|
||||
p % coord % universe = BASE_UNIVERSE
|
||||
|
||||
! loop over horizontal rays
|
||||
do while(xyz(2) > last_y_coord)
|
||||
do y=1, img%height
|
||||
do x=1, img%width
|
||||
|
||||
! initialize the particle and set starting coordinate and direction
|
||||
call initialize_particle()
|
||||
call deallocate_coord(p % coord0 % next)
|
||||
p % coord => p % coord0
|
||||
|
||||
p % coord % xyz = xyz
|
||||
p % coord % uvw = (/ 1, 0, 0 /)
|
||||
call find_cell(found_cell)
|
||||
|
||||
! write starting coordinate to file
|
||||
write(UNIT=UNIT_PLOT) p % coord % xyz
|
||||
|
||||
! Find cell that particle is currently in
|
||||
call find_cell(found_cell)
|
||||
|
||||
! =======================================================================
|
||||
! MOVE PARTICLE FORWARD TO NEXT CELL
|
||||
|
||||
if (.not. found_cell) then
|
||||
! Clear any coordinates beyond first level
|
||||
call deallocate_coord(p % coord0 % next)
|
||||
p % coord => p % coord0
|
||||
|
||||
distance = INFINITY
|
||||
univ => universes(BASE_UNIVERSE)
|
||||
do i = 1, univ % n_cells
|
||||
p % coord0 % xyz = xyz
|
||||
p % coord0 % cell = univ % cells(i)
|
||||
|
||||
call distance_to_boundary(d, surface_crossed, lattice_crossed)
|
||||
if (d < distance) then
|
||||
! Check to make sure particle is actually going into this cell
|
||||
! by moving it slightly forward and seeing if the cell contains
|
||||
! that coordinate
|
||||
|
||||
p % coord0 % xyz = p % coord0 % xyz + (d + TINY_BIT) * p % coord0 % uvw
|
||||
|
||||
c => cells(p % coord0 % cell)
|
||||
if (.not. cell_contains(c)) cycle
|
||||
|
||||
! Set new distance and retain pointer to this cell
|
||||
distance = d
|
||||
enter_surface = surface_crossed
|
||||
end if
|
||||
end do
|
||||
|
||||
! No cell was found on this horizontal ray
|
||||
if (distance == INFINITY) then
|
||||
p % coord0 % xyz(1) = last_x_coord
|
||||
write(UNIT_PLOT) p % coord0 % xyz, 0
|
||||
|
||||
! Move to next horizontal ray
|
||||
xyz(2) = xyz(2) - pixel
|
||||
cycle
|
||||
end if
|
||||
|
||||
! Write coordinate where next cell begins
|
||||
p % coord0 % xyz = xyz + distance * p % coord0 % uvw
|
||||
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
|
||||
|
||||
! Process surface crossing for next cell
|
||||
p % coord0 % cell = NONE
|
||||
p % surface = -enter_surface
|
||||
call cross_surface(enter_surface)
|
||||
end if
|
||||
|
||||
! =======================================================================
|
||||
! 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(distance, surface_crossed, lattice_crossed)
|
||||
|
||||
! Advance particle
|
||||
coord => p % coord0
|
||||
do while (associated(coord))
|
||||
coord % xyz = coord % xyz + distance * coord % uvw
|
||||
coord => coord % next
|
||||
end do
|
||||
|
||||
! If next boundary crossing is out of range of the plot, only include
|
||||
! the visible portion and move to next horizontal ray
|
||||
if (p % coord0 % xyz(1) >= last_x_coord) then
|
||||
p % alive = .false.
|
||||
p % coord0 % xyz(1) = last_x_coord
|
||||
|
||||
! If there is no cell beyond this boundary, mark it as cell 0
|
||||
if (distance == INFINITY) p % coord % cell = 0
|
||||
|
||||
! Write ending coordinates to file
|
||||
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, last_cell
|
||||
|
||||
p % coord % cell = 0
|
||||
if (lattice_crossed /= NONE) then
|
||||
p % surface = NONE
|
||||
call cross_lattice(lattice_crossed)
|
||||
if (.not. found_cell) then
|
||||
r = 0
|
||||
g = 0
|
||||
b = 0
|
||||
else
|
||||
c => cells(p % coord % cell)
|
||||
if (pl % color == PLOT_COLOR_MATS) then
|
||||
r = materials(c % material) % rgb(1)
|
||||
g = materials(c % material) % rgb(2)
|
||||
b = materials(c % material) % rgb(3)
|
||||
else if (pl % color == PLOT_COLOR_CELLS) then
|
||||
r = c % rgb(1)
|
||||
g = c % rgb(2)
|
||||
b = c % rgb(3)
|
||||
else
|
||||
p % surface = surface_crossed
|
||||
call cross_surface(last_cell)
|
||||
|
||||
! Since boundary conditions are disabled in plotting mode, we need
|
||||
! to manually add the last segment
|
||||
if (surfaces(abs(surface_crossed)) % bc /= BC_TRANSMIT) then
|
||||
p % coord0 % xyz(1) = last_x_coord
|
||||
write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
|
||||
exit
|
||||
end if
|
||||
r = 0
|
||||
g = 0
|
||||
b = 0
|
||||
end if
|
||||
end if
|
||||
|
||||
call set_pixel(img, x, y, r, g, b)
|
||||
|
||||
p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
|
||||
end do
|
||||
|
||||
p % coord0 % xyz(in_i) = xyz(in_i)
|
||||
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel
|
||||
end do
|
||||
|
||||
call output_ppm(pl,img)
|
||||
|
||||
call deallocate_image(img)
|
||||
|
||||
end subroutine create_ppm
|
||||
|
||||
|
||||
!===============================================================================
|
||||
! output_ppm
|
||||
!===============================================================================
|
||||
subroutine output_ppm(pl,img)
|
||||
|
||||
type(Plot), pointer :: pl
|
||||
type(Image), intent(in) :: img
|
||||
|
||||
integer :: i, j
|
||||
character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
|
||||
|
||||
path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm"
|
||||
open(UNIT=UNIT_PLOT, FILE=path_plot)
|
||||
|
||||
write(UNIT_PLOT, '(A2)') 'P6'
|
||||
write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height
|
||||
write(UNIT_PLOT, '(A)') '255'
|
||||
|
||||
do j=1, img%height
|
||||
do i=1, img%width
|
||||
write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), &
|
||||
achar(img%green(i,j)), &
|
||||
achar(img%blue(i,j))
|
||||
end do
|
||||
|
||||
! Move y-coordinate to next position
|
||||
xyz(2) = xyz(2) - pixel
|
||||
end do
|
||||
|
||||
! Close plot file
|
||||
close(UNIT=UNIT_PLOT)
|
||||
|
||||
end subroutine output_ppm
|
||||
|
||||
end subroutine run_plot
|
||||
|
||||
end module plot
|
||||
end module plotter
|
||||
|
|
|
|||
31
src/plot_header.F90
Normal file
31
src/plot_header.F90
Normal file
|
|
@ -0,0 +1,31 @@
|
|||
module plot_header
|
||||
|
||||
implicit none
|
||||
|
||||
!===============================================================================
|
||||
! PLOT hold plot information
|
||||
!===============================================================================
|
||||
|
||||
type Plot
|
||||
integer :: id ! Unique ID
|
||||
integer :: type ! Type
|
||||
integer :: color ! quantity to color regions by
|
||||
real(8) :: origin(3) ! xyz center of plot location
|
||||
real(8) :: aspect ! spacing between rays in raytracer
|
||||
real(8) :: width(3) ! xyz widths of plot
|
||||
integer :: basis ! direction of plot slice
|
||||
integer :: pixels(2) ! pixel width/height of plot slice
|
||||
end type Plot
|
||||
|
||||
integer :: PLOT_TYPE_SLICE = 1
|
||||
integer :: PLOT_TYPE_POINTS = 2
|
||||
|
||||
integer :: PLOT_BASIS_XY = 1
|
||||
integer :: PLOT_BASIS_XZ = 2
|
||||
integer :: PLOT_BASIS_YZ = 3
|
||||
|
||||
integer :: PLOT_COLOR_CELLS = 1
|
||||
integer :: PLOT_COLOR_MATS = 2
|
||||
|
||||
|
||||
end module plot_header
|
||||
101
src/ppmlib.F90
Normal file
101
src/ppmlib.F90
Normal file
|
|
@ -0,0 +1,101 @@
|
|||
module ppmlib
|
||||
|
||||
implicit none
|
||||
|
||||
!===============================================================================
|
||||
! Image holds RGB information for output PPM image
|
||||
!===============================================================================
|
||||
|
||||
type Image
|
||||
integer, dimension(:,:), pointer :: red, green, blue
|
||||
integer :: width, height
|
||||
end type Image
|
||||
|
||||
contains
|
||||
|
||||
subroutine init_image(img)
|
||||
|
||||
type(Image), intent(out) :: img
|
||||
|
||||
nullify(img % red)
|
||||
nullify(img % green)
|
||||
nullify(img % blue)
|
||||
|
||||
img % width = 0
|
||||
img % height = 0
|
||||
|
||||
end subroutine init_image
|
||||
|
||||
subroutine allocate_image(img, w, h)
|
||||
|
||||
type(Image), intent(inout) :: img
|
||||
integer, intent(in) :: w, h
|
||||
|
||||
allocate(img % red(w, h))
|
||||
allocate(img % green(w, h))
|
||||
allocate(img % blue(w, h))
|
||||
|
||||
img % width = w
|
||||
img % height = h
|
||||
|
||||
end subroutine allocate_image
|
||||
|
||||
subroutine deallocate_image(img)
|
||||
|
||||
type(Image) :: img
|
||||
|
||||
if ( associated(img % red) ) deallocate(img % red)
|
||||
if ( associated(img % green) ) deallocate(img % green)
|
||||
if ( associated(img % blue) ) deallocate(img % blue)
|
||||
end subroutine deallocate_image
|
||||
|
||||
|
||||
function inside_image(img, x, y) result(inside)
|
||||
|
||||
type(Image), intent(in) :: img
|
||||
integer, intent(in) :: x, y
|
||||
logical :: inside
|
||||
|
||||
inside = .false.
|
||||
|
||||
if ( ( x < img % width) .and. &
|
||||
( y < img % height ) .and. &
|
||||
( x >= 0 ) .and. &
|
||||
( y >= 0 ) ) inside = .true.
|
||||
|
||||
end function inside_image
|
||||
|
||||
|
||||
function valid_image(img) result(valid)
|
||||
|
||||
type(Image), intent(in) :: img
|
||||
logical :: valid
|
||||
|
||||
valid = .false.
|
||||
|
||||
if ( img % width == 0 ) return
|
||||
if ( img % height == 0 ) return
|
||||
if ( .not. associated(img % red) .or. &
|
||||
.not. associated(img % green) .or. &
|
||||
.not. associated(img % blue) ) return
|
||||
|
||||
valid = .true.
|
||||
|
||||
end function valid_image
|
||||
|
||||
|
||||
subroutine set_pixel(img, x, y, r, g, b)
|
||||
|
||||
type(Image), intent(inout) :: img
|
||||
integer, intent(in) :: x, y
|
||||
integer, intent(in) :: r, g, b
|
||||
|
||||
if ( inside_image(img, x, y) .and. valid_image(img)) then
|
||||
img % red(x+1,y+1) = mod(abs(r), 256)
|
||||
img % green(x+1, y+1) = mod(abs(g), 256)
|
||||
img % blue(x+1, y+1) = mod(abs(b), 256)
|
||||
end if
|
||||
end subroutine set_pixel
|
||||
|
||||
|
||||
end module ppmlib
|
||||
|
|
@ -3,9 +3,17 @@
|
|||
|
||||
<options rootname="plot" />
|
||||
|
||||
<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" default="'xy'" />
|
||||
<variable name="pixel_" tag="pixel" type="double" default="0.01" />
|
||||
<typedef name="plot_xml">
|
||||
<component name="id" type="integer" />
|
||||
<component name="type" type="word" length="10" default="''"/>
|
||||
<component name="color" type="word" length="10" default="'cell'"/>
|
||||
<component name="origin" type="double-array" />
|
||||
<component name="width" type="double-array" />
|
||||
<component name="basis" type="word" length="3" default="'xy'" />
|
||||
<component name="pixels" type="integer-array"/>
|
||||
<component name="aspect" type="double" default="0.1" />
|
||||
</typedef>
|
||||
|
||||
<variable name="plot_" tag="plot" type="plot_xml" dimension="1" />
|
||||
|
||||
</template>
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue