diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index c184ea2134..c42098a765 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -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
diff --git a/src/finalize.F90 b/src/finalize.F90
index d601c2e592..9fd4bcb0af 100644
--- a/src/finalize.F90
+++ b/src/finalize.F90
@@ -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()
diff --git a/src/global.F90 b/src/global.F90
index 83e31758e7..3dd62515fe 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -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
diff --git a/src/initialize.F90 b/src/initialize.F90
index a1518cd218..fc7569baf1 100644
--- a/src/initialize.F90
+++ b/src/initialize.F90
@@ -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()
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 13ff553a92..8eb5e29fc1 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -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)
diff --git a/src/main.F90 b/src/main.F90
index 5e5ea89ba1..e7456f5055 100644
--- a/src/main.F90
+++ b/src/main.F90
@@ -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
diff --git a/src/material_header.F90 b/src/material_header.F90
index cbfd917498..e8ccfff344 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -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
diff --git a/src/mesh.F90 b/src/mesh.F90
index f24360c446..e69f6ab5cb 100644
--- a/src/mesh.F90
+++ b/src/mesh.F90
@@ -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
diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90
index 96cc587a6e..9cb371872f 100644
--- a/src/particle_restart_write.F90
+++ b/src/particle_restart_write.F90
@@ -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', &
diff --git a/src/physics.F90 b/src/physics.F90
index aa2f881bdc..cee549aa71 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -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
diff --git a/src/plot.F90 b/src/plot.F90
index 6252b6c6b9..f5b20c8515 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -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
diff --git a/src/plot_header.F90 b/src/plot_header.F90
index 6956bfaf93..8030f61d68 100644
--- a/src/plot_header.F90
+++ b/src/plot_header.F90
@@ -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
diff --git a/src/templates/plots_t.xml b/src/templates/plots_t.xml
index b3b86fc44e..8781abf016 100644
--- a/src/templates/plots_t.xml
+++ b/src/templates/plots_t.xml
@@ -25,6 +25,7 @@
+
diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml
index 3afa61745e..2314c7c3d8 100644
--- a/src/templates/settings_t.xml
+++ b/src/templates/settings_t.xml
@@ -54,6 +54,7 @@
+