Initial add of support for serpent-style plots

This commit is contained in:
nhorelik 2013-05-06 09:37:53 -04:00
parent e950c52067
commit 29fc61bb2a
14 changed files with 240 additions and 24 deletions

View file

@ -128,6 +128,7 @@ finalize.o: cmfd_output.o
finalize.o: global.o
finalize.o: hdf5_interface.o
finalize.o: output.o
finalize.o: plot.o
finalize.o: tally.o
fission.o: ace_header.o
@ -297,6 +298,7 @@ physics.o: mesh.o
physics.o: output.o
physics.o: particle_header.o
physics.o: particle_restart_write.o
physics.o: plot.o
physics.o: random_lcg.o
physics.o: search.o
physics.o: string.o
@ -307,15 +309,16 @@ plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
plot.o: mesh.o
plot.o: output.o
plot.o: particle_header.o
plot.o: plot_header.o
plot.o: ppmlib.o
plot.o: progress_header.o
plot.o: source.o
plot.o: string.o
plot_header.o: constants.o
plot_header.o: mesh_header.o
random_lcg.o: global.o

View file

@ -5,6 +5,7 @@ module finalize
# endif
use global
use output, only: print_runtime, print_results, write_tallies
use plot, only: finalize_rxn_plots
use tally, only: tally_statistics
#ifdef MPI
@ -52,6 +53,9 @@ contains
call print_results()
end if
! Finish reaction rate plots
if (master .and. rxn_plots) call finalize_rxn_plots()
! deallocate arrays
call free_memory()

View file

@ -396,14 +396,23 @@ contains
if (allocated(materials)) deallocate(materials)
! Deallocate plots
if (allocated(plots)) then
do i = 1, size(plots)
if allocated(plots(i) % colors) deallocate(plots(i) % colors)
if allocated(plots(i) % rxncnt) deallocate(plots(i) % rxncnt)
if (allocated(plots(i) % colors)) deallocate(plots(i) % colors)
if (allocated(plots(i) % fisswgt)) deallocate(plots(i) % fisswgt)
if (allocated(plots(i) % fluxwgt)) deallocate(plots(i) % fluxwgt)
if (allocated(plots(i) % pixmesh % dimension)) &
deallocate(plots(i) % pixmesh % dimension)
if (allocated(plots(i) % pixmesh % lower_left)) &
deallocate(plots(i) % pixmesh % lower_left)
if (allocated(plots(i) % pixmesh % upper_right)) &
deallocate(plots(i) % pixmesh % upper_right)
if (allocated(plots(i) % pixmesh % width)) &
deallocate(plots(i) % pixmesh % width)
end do
deallocate(plots)
end if
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) then
! First call the clear routines

View file

@ -76,7 +76,7 @@ contains
! Read plots.xml if it exists -- this has to be done separate from the other
! XML files because we need the PRNG to be initialized first
if (run_mode == MODE_PLOTTING) call read_plots_xml()
if (run_mode == MODE_PLOTTING .or. rxn_plots) call read_plots_xml()
! Set up universe structures
call prepare_universes()
@ -98,6 +98,11 @@ contains
call read_xs()
call time_read_xs % stop()
TODO
! Set material fissionable flag
matnuc => nuclides(mat % nuclide(j))
if (matnuc % fissionable) mat % fissionable = .true.
! Construct unionized energy grid from cross-sections
if (grid_method == GRID_UNION) then
call time_unionize % start()

View file

@ -570,6 +570,13 @@ contains
#endif
end if
! Check flag to enable rxn plot scoring
call lower_case(rxn_plots_)
if (rxn_plots_ == 'true' .or. rxn_plots == '1') then
rxn_plots = .true.
end if
end subroutine read_settings_xml
!===============================================================================
@ -1011,6 +1018,7 @@ contains
type(ListChar) :: list_names ! temporary list of nuclide names
type(ListReal) :: list_density ! temporary list of nuclide densities
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: matnuc => null()
type(nuclide_xml), pointer :: nuc => null()
type(sab_xml), pointer :: sab => null()
@ -1257,6 +1265,7 @@ contains
! Copy name and atom/weight percent
mat % names(j) = name
mat % atom_density(j) = list_density % get_item(j)
end do ALL_NUCLIDES
! Check to make sure either all atom percents or all weight percents are
@ -2184,6 +2193,8 @@ contains
pl % type = PLOT_TYPE_SLICE
case ("voxel")
pl % type = PLOT_TYPE_VOXEL
case ("rxn")
pl % type = PLOT_TYPE_RXNRATE
case default
message = "Unsupported plot type '" // trim(plot_(i) % type) &
// "' in plot " // trim(to_str(pl % id))
@ -2192,7 +2203,7 @@ contains
! Set output file path
select case (pl % type)
case (PLOT_TYPE_SLICE)
case (PLOT_TYPE_SLICE, PLOT_TYPE_RXNRATE)
pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
"_" // trim(plot_(i) % filename) // ".ppm"
case (PLOT_TYPE_VOXEL)
@ -2201,7 +2212,8 @@ contains
end select
! Copy plot pixel size
if (pl % type == PLOT_TYPE_SLICE) then
if (pl % type == PLOT_TYPE_SLICE .or. &
pl % type == PLOT_TYPE_RXNRATE) then
if (size(plot_(i) % pixels) == 2) then
pl % pixels(1) = plot_(i) % pixels(1)
pl % pixels(2) = plot_(i) % pixels(2)
@ -2239,7 +2251,8 @@ contains
end if
! Copy plot basis
if (pl % type == PLOT_TYPE_SLICE) then
if (pl % type == PLOT_TYPE_SLICE .or. &
pl % type == PLOT_TYPE_RXNRATE) then
select case (plot_(i) % basis)
case ("xy")
pl % basis = PLOT_BASIS_XY
@ -2264,7 +2277,8 @@ contains
end if
! Copy plotting width
if (pl % type == PLOT_TYPE_SLICE) then
if (pl % type == PLOT_TYPE_SLICE .or. &
pl % type == PLOT_TYPE_RXNRATE) then
if (size(plot_(i) % width) == 2) then
pl % width(1) = plot_(i) % width(1)
pl % width(2) = plot_(i) % width(2)
@ -2412,6 +2426,48 @@ contains
end if
! Initialize attributes specific to reaction rate plots
if (pl % type == PLOT_TYPE_RXNRATE) then
! Copy rxn rate type
select case (plot_(i) % rrtype)
case ("fission")
pl % rrtype = PLOT_RXN_FISSION
case ("absorption")
pl % rrtype = PLOT_RXN_ABSORPTION
case ("fluxthermal")
pl % rrtype = PLOT_RXN_FLUX_THERMAL
case ("fluxfast")
pl % rrtype = PLOT_RXN_FLUX_FAST
case default
message = "Unsupported plot rxn type '" // trim(plot_(i) % rrtype) &
// "' in plot " // trim(to_str(pl % id))
call fatal_error()
end select
! Allocate and set pixel mesh parameters
allocate(pl % pixmesh % dimension(2))
allocate(pl % pixmesh % lower_left(2))
allocate(pl % pixmesh % upper_right(2))
allocate(pl % pixmesh % width(2))
pl % pixmesh % n_dimension = 2
pl % pixmesh % dimension(1) = pl % pixels(1)
pl % pixmesh % dimension(2) = pl % pixels(2)
pl % pixmesh % lower_left(1) = pl % origin(1) - pl % width(1) / 2.0
pl % pixmesh % lower_left(2) = pl % origin(2) - pl % width(2) / 2.0
pl % pixmesh % upper_right(1) = pl % origin(1) + pl % width(1) / 2.0
pl % pixmesh % upper_right(2) = pl % origin(2) + pl % width(2) / 2.0
pl % pixmesh % width(1) = pl % width(1) / dble(pl % pixels(1))
pl % pixmesh % width(2) = pl % width(2) / dble(pl % pixels(2))
! Initialize reaction counter
allocate(pl % fisswgt(pl % pixels(1) * pl % pixels(2)))
allocate(pl % fluxwgt(pl % pixels(1) * pl % pixels(2)))
pl % fisswgt = ZERO
pl % fluxwgt = ZERO
end if
! Add plot to dictionary
call plot_dict % add_key(pl % id, i)

View file

@ -7,7 +7,7 @@ program main
use global
use initialize, only: initialize_run
use particle_restart, only: run_particle_restart
use plot, only: run_plot
use plot, only: run_plots
implicit none
@ -21,7 +21,7 @@ program main
case (MODE_EIGENVALUE)
call run_eigenvalue()
case (MODE_PLOTTING)
call run_plot()
call run_plots()
case (MODE_TALLIES)
! For tallies-only mode, we just skip straight to finalize_run to write out
! the tally results

View file

@ -21,6 +21,9 @@ module material_header
! Temporary names read during initialization
character(12), allocatable :: names(:) ! isotope names
character(12), allocatable :: sab_names(:) ! name of S(a,b) table
logical :: fissionable = .false.
end type Material
end module material_header

View file

@ -2,7 +2,7 @@ module mesh
use constants
use global
use mesh_header
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
use search, only: binary_search

View file

@ -90,7 +90,7 @@ contains
type(Bank), pointer :: src => null()
! set up file name
filename = path_output // 'particle_'//trim(to_str(rank))//'.binary'
filename = trim(path_output) // 'particle_'//trim(to_str(rank))//'.binary'
! create hdf5 file
open(UNIT=UNIT_PARTICLE, FILE=filename, STATUS='replace', &

View file

@ -14,6 +14,7 @@ module physics
use material_header, only: Material
use mesh, only: get_mesh_indices
use output, only: write_message
use plot, only: score_rxn_rate_plots
use particle_header, only: LocalCoord
use particle_restart_write, only: write_particle_restart
use random_lcg, only: prn
@ -140,6 +141,8 @@ contains
p % surface = NONE
call collision()
if (rxn_plots .and. active_batches) call score_rxn_rate_plots()
! Save coordinates for tallying purposes
p % last_xyz = p % coord0 % xyz

View file

@ -5,13 +5,13 @@ module plot
use geometry, only: find_cell
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use mesh, only: get_mesh_bin, mesh_indices_to_bin
use output, only: write_message
use particle_header, only: deallocate_coord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
use progress_header
use source, only: initialize_particle
use string, only: to_str
implicit none
@ -19,10 +19,10 @@ module plot
contains
!===============================================================================
! RUN_PLOT controls the logic for making one or many plots
! RUN_PLOTs controls the logic for making one or many plots in MODE_PLOTTING
!===============================================================================
subroutine run_plot()
subroutine run_plots()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl => null()
@ -43,22 +43,96 @@ contains
end if
end do
end subroutine run_plot
end subroutine run_plots
!===============================================================================
! SCORE_RXN_RATE_PLOTS
!===============================================================================
subroutine score_rxn_rate_plots()
integer :: i ! loop index for plots
integer :: n
integer :: bin
type(ObjectPlot), pointer :: pl => null()
type(StructuredMesh), pointer :: m => null()
do i = 1, n_plots
pl => plots(i)
if (.not. pl % type == PLOT_TYPE_RXNRATE) cycle
if (.not. nuclides(p % event_nuclide) % fissionable) then
select case(pl % rrtype)
case (PLOT_RXN_FLUX_THERMAL)
if (.not. p % E <= 0.625e-6) cycle
case (PLOT_RXN_FLUX_FAST)
if (.not. p % E >= 0.625e-6) cycle
end select
end if
m => pl % pixmesh
call get_mesh_bin(m, p % coord0 % xyz, bin)
if (bin == NO_BIN_FOUND) cycle
if (nuclides(p % event_nuclide) % fissionable) then
! pl % fisswgt(bin) = pl % fisswgt(bin) + keff * p % wgt_bank
n = nuclides(p % event_nuclide) % index_fission(1)
pl % fisswgt(bin) = pl% fisswgt(bin) + p % last_wgt * &
nuclides(p % event_nuclide) % reactions(n) % Q_value
else
pl % fluxwgt(bin) = pl % fluxwgt(bin) + p % last_wgt / material_xs % total
end if
end do
end subroutine score_rxn_rate_plots
!===============================================================================
! FINALIZE_RXN_PLOTS
!===============================================================================
subroutine finalize_rxn_plots()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl => null()
do i = 1, n_plots
pl => plots(i)
if (.not. pl % type == PLOT_TYPE_RXNRATE) cycle
pl % fisswgt = pl % fisswgt / maxval(pl % fisswgt)
pl % fluxwgt = pl % fluxwgt / maxval(pl % fluxwgt)
! Display output message
message = "Processing rxn plot " // trim(to_str(pl % id)) // "..."
call write_message(5)
call create_ppm(pl, rxn_overlay=.true.)
end do
end subroutine finalize_rxn_plots
!===============================================================================
! POSITION_RGB computes the red/green/blue values for a given plot with the
! current particle's position
!===============================================================================
subroutine position_rgb(pl, rgb, id)
subroutine position_rgb(pl, rgb, id, fissionable)
type(ObjectPlot), pointer, intent(in) :: pl
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
logical, intent(out), optional :: fissionable
logical :: found_cell
type(Cell), pointer :: c => null()
if (present(fissionable)) fissionable = .false.
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
@ -69,9 +143,9 @@ contains
rgb = pl % not_found % rgb
id = -1
else
c => cells(p % coord % cell)
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord % cell)
id = materials(c % material) % id
if (c % material == MATERIAL_VOID) then
! By default, color void cells white
@ -87,6 +161,9 @@ contains
rgb = 0
id = -1
end if
if (present(fissionable)) then
fissionable = materials(c % material) % fissionable
end if
end if
end subroutine position_rgb
@ -96,20 +173,32 @@ contains
! specification in the portable pixmap format (PPM)
!===============================================================================
subroutine create_ppm(pl)
subroutine create_ppm(pl, rxn_overlay)
type(ObjectPlot), pointer :: pl
logical, optional :: rxn_overlay
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
integer :: ijk(2) ! pixmesh indices
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id
integer :: bin ! pixmesh results bin
real(8) :: rxnval ! pixmesh results value
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
logical :: rxn = .false.
logical :: fissionable
type(Image) :: img
type(ProgressBar) :: progress
type(StructuredMesh), pointer :: m => null()
if (present(rxn_overlay) .and. rxn_overlay) then
rxn = .true.
m => pl % pixmesh
end if
! Initialize and allocate space for image
call init_image(img)
@ -142,8 +231,13 @@ contains
end if
! allocate and initialize particle
! We can't depend on source.F90 for initialize_particle due to
! circular dependencies (since physics.F90 depends on this file)
! To fix this, transport and physics could be
allocate(p)
call initialize_particle()
allocate(p % coord0)
p % coord0 % universe = BASE_UNIVERSE
p % coord => p % coord0
p % coord % xyz = xyz
p % coord % uvw = (/ 0.5, 0.5, 0.5 /)
p % coord % universe = BASE_UNIVERSE
@ -152,8 +246,26 @@ contains
call progress % set_value(dble(y)/dble(img % height)*100.)
do x = 1, img % width
! get pixel color
call position_rgb(pl, rgb, id)
! Get pixel color
call position_rgb(pl, rgb, id, fissionable)
! Add rxn rate overlay if appropriate
if (rxn) then
ijk(1) = x
ijk(2) = y
bin = mesh_indices_to_bin(m, ijk)
if (fissionable) then
rxnval = pl % fisswgt(bin)
rgb(1) = min(int(1010.*rxnval),255)
rgb(2) = int(255.*rxnval)
rgb(3) = int(255.*rxnval)
else
rxnval = pl % fluxwgt(bin)
rgb(1) = int(255.*rxnval)
rgb(2) = int(255.*rxnval)
rgb(3) = min(int(1010.*rxnval),255)
end if
end if
! Create a pixel at (x,y) with color (r,g,b)
call set_pixel(img, x, y, rgb(1), rgb(2), rgb(3))
@ -238,8 +350,13 @@ contains
ll = pl % origin - pl % width / 2.0
! allocate and initialize particle
! We can't depend on source.F90 for initialize_particle due to
! circular dependencies (since physics.F90 depends on this file)
! To fix this, transport and physics could be
allocate(p)
call initialize_particle()
allocate(p % coord0)
p % coord0 % universe = BASE_UNIVERSE
p % coord => p % coord0
p % coord0 % xyz = ll
p % coord0 % uvw = (/ 0.5, 0.5, 0.5 /)
p % coord0 % universe = BASE_UNIVERSE

View file

@ -1,6 +1,7 @@
module plot_header
use constants
use mesh_header, only: StructuredMesh
implicit none
@ -27,11 +28,16 @@ module plot_header
integer :: pixels(3) ! pixel width/height of plot slice
type(ObjectColor) :: not_found ! color for positions where no cell found
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
type(StructuredMesh) :: pixmesh ! pixmesh for reaction rate plots
real(8), allocatable :: fisswgt(:)! Weights for rxn plot fissionable pixels
real(8), allocatable :: fluxwgt(:)! Weights for rxn plot non-fiss pix
integer :: rrtype ! particle event type to score
end type ObjectPlot
! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
integer, parameter :: PLOT_TYPE_VOXEL = 2
integer, parameter :: PLOT_TYPE_RXNRATE = 3
! Plot basis plane
integer, parameter :: PLOT_BASIS_XY = 1
@ -42,4 +48,12 @@ module plot_header
integer, parameter :: PLOT_COLOR_CELLS = 1
integer, parameter :: PLOT_COLOR_MATS = 2
! Reaction rate plot types
integer, parameter :: PLOT_RXN_FLUX_THERMAL = 1
integer, parameter :: PLOT_RXN_FLUX_FAST = 2
integer, parameter :: PLOT_RXN_FISSION = 3
integer, parameter :: PLOT_RXN_ABSORPTION = 4
contains
end module plot_header

View file

@ -25,6 +25,7 @@
<component name="background" type="integer-array"/>
<component name="col_spec_" tag="col_spec" type="col_spec_xml" dimension="1" />
<component name="mask_" tag="mask" type="mask_xml" dimension="1" />
<component name="rrtype" type="word" length="15" default="'fission'" />
</typedef>
<variable name="plot_" tag="plot" type="plot_xml" dimension="1" />

View file

@ -54,6 +54,7 @@
<variable name="output_path_" tag="output_path" type="word" length="255" />
<variable name="ptables_" tag="ptables" type="word" length="5" />
<variable name="run_cmfd_" tag="run_cmfd" type="word" length="5" />
<variable name="rxn_plots_" tag="enable_rxn_plots" type="word" length="5" />
<variable name="seed_" tag="seed" type="integer" />
<variable name="source_" tag="source" type="source_xml" />
<variable name="state_point_" tag="state_point" type="statepoint_xml" dimension="1" />