From 29fc61bb2ae555ac5d3feae600117985cffddc06 Mon Sep 17 00:00:00 2001 From: nhorelik Date: Mon, 6 May 2013 09:37:53 -0400 Subject: [PATCH] Initial add of support for serpent-style plots --- src/DEPENDENCIES | 5 +- src/finalize.F90 | 4 + src/global.F90 | 15 +++- src/initialize.F90 | 7 +- src/input_xml.F90 | 64 ++++++++++++++- src/main.F90 | 4 +- src/material_header.F90 | 3 + src/mesh.F90 | 2 +- src/particle_restart_write.F90 | 2 +- src/physics.F90 | 3 + src/plot.F90 | 139 ++++++++++++++++++++++++++++++--- src/plot_header.F90 | 14 ++++ src/templates/plots_t.xml | 1 + src/templates/settings_t.xml | 1 + 14 files changed, 240 insertions(+), 24 deletions(-) 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 @@ +