From 40a931f97bfdc03380820f86cf9243f218baf99b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 20 Feb 2017 13:17:50 -0600 Subject: [PATCH] Parallelize slice plots --- CMakeLists.txt | 1 - src/plot.F90 | 118 +++++++++++++++++++++-------------------- src/ppmlib.F90 | 127 --------------------------------------------- src/ppmlib.LICENSE | 4 -- 4 files changed, 58 insertions(+), 192 deletions(-) delete mode 100644 src/ppmlib.F90 delete mode 100644 src/ppmlib.LICENSE diff --git a/CMakeLists.txt b/CMakeLists.txt index 098de3de2a..3ecdc1e9bc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -318,7 +318,6 @@ set(LIBOPENMC_FORTRAN_SRC src/physics_mg.F90 src/plot.F90 src/plot_header.F90 - src/ppmlib.F90 src/product_header.F90 src/progress_header.F90 src/random_lcg.F90 diff --git a/src/plot.F90 b/src/plot.F90 index dd0d4bfbd8..86d303b3dd 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -11,14 +11,19 @@ module plot use output, only: write_message, time_stamp use particle_header, only: LocalCoord, Particle use plot_header - use ppmlib, only: Image, init_image, allocate_image, & - deallocate_image, set_pixel use progress_header, only: ProgressBar use string, only: to_str use hdf5 implicit none + private + + public :: run_plot + + integer, parameter :: RED = 1 + integer, parameter :: GREEN = 2 + integer, parameter :: BLUE = 3 contains @@ -41,7 +46,7 @@ contains call create_ppm(pl) else if (pl % type == PLOT_TYPE_VOXEL) then ! create dump for 3D silomesh utility script - call create_3d_dump(pl) + call create_voxel(pl) end if end associate end do @@ -117,19 +122,22 @@ contains integer :: x, y ! pixel location integer :: rgb(3) ! colors (red, green, blue) from 0-255 integer :: id + integer :: height, width real(8) :: in_pixel real(8) :: out_pixel real(8) :: xyz(3) - type(Image) :: img + integer, allocatable :: data(:,:,:) type(Particle) :: p - type(ProgressBar) :: progress - ! Initialize and allocate space for image - call init_image(img) - call allocate_image(img, pl % pixels(1), pl % pixels(2)) + width = pl % pixels(1) + height = pl % pixels(2) - in_pixel = pl % width(1)/dble(pl % pixels(1)) - out_pixel = pl % width(2)/dble(pl % pixels(2)) + in_pixel = pl % width(1)/dble(width) + out_pixel = pl % width(2)/dble(height) + + ! Allocate and initialize results array + allocate(data(3, width, height)) + data(:,:,:) = 0 if (pl % basis == PLOT_BASIS_XY) then in_i = 1 @@ -157,36 +165,28 @@ contains p % coord(1) % uvw = [ HALF, HALF, HALF ] p % coord(1) % universe = BASE_UNIVERSE - do y = 1, img % height - call progress % set_value(dble(y)/dble(img % height)*100) - do x = 1, img % width +!$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data) + do y = 1, height + ! Set y coordinate + p % coord(1) % xyz(out_i) = xyz(out_i) - out_pixel*(y - 1) + do x = 1, width + ! Set x coordinate + p % coord(1) % xyz(in_i) = xyz(in_i) + in_pixel*(x - 1) ! get pixel color call position_rgb(p, pl, rgb, id) ! Create a pixel at (x,y) with color (r,g,b) - call set_pixel(img, x-1, y-1, rgb(1), rgb(2), rgb(3)) - - ! Advance pixel in first direction - p % coord(1) % xyz(in_i) = p % coord(1) % xyz(in_i) + in_pixel + data(:, x, y) = rgb end do - - ! Advance pixel in second direction - p % coord(1) % xyz(in_i) = xyz(in_i) - p % coord(1) % xyz(out_i) = p % coord(1) % xyz(out_i) - out_pixel end do +!$omp end parallel do ! Draw tally mesh boundaries on the image if requested - if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, img) + if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, data) ! Write out the ppm to a file - call output_ppm(pl, img) - - ! Free up space - call deallocate_image(img) - - ! Clear particle - call p % clear() + call output_ppm(pl, data) end subroutine create_ppm @@ -194,13 +194,13 @@ contains ! DRAW_MESH_LINES draws mesh line boundaries on an image !=============================================================================== - subroutine draw_mesh_lines(pl, img) + subroutine draw_mesh_lines(pl, data) type(ObjectPlot), intent(in) :: pl - type(Image), intent(inout) :: img + integer, intent(inout) :: data(:,:,:) logical :: in_mesh integer :: out_, in_ ! pixel location - integer :: r, g, b ! RGB color for meshlines pixels + integer :: rgb(3) ! RGB color for meshlines pixels integer :: outrange(2), inrange(2) ! range of pixel locations integer :: i, j ! loop indices integer :: plus @@ -214,9 +214,7 @@ contains real(8) :: xyz_ll(3) ! lower left xyz real(8) :: xyz_ur(3) ! upper right xyz - r = pl % meshlines_color % rgb(1) - g = pl % meshlines_color % rgb(2) - b = pl % meshlines_color % rgb(3) + rgb(:) = pl % meshlines_color % rgb select case (pl % basis) case(PLOT_BASIS_XY) @@ -260,30 +258,30 @@ contains ! map the xyz ranges to pixel ranges frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer) - outrange(1) = int(frac * real(img % width, 8)) + outrange(1) = int(frac * real(pl % pixels(1), 8)) frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) - outrange(2) = int(frac * real(img % width, 8)) + outrange(2) = int(frac * real(pl % pixels(1), 8)) frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(1) = int((ONE - frac) * real(img % height, 8)) + inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8)) frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(2) = int((ONE - frac) * real(img % height, 8)) + inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8)) ! draw lines do out_ = outrange(1), outrange(2) do plus = 0, pl % meshlines_width - call set_pixel(img, out_, inrange(1) + plus, r, g, b) - call set_pixel(img, out_, inrange(2) + plus, r, g, b) - call set_pixel(img, out_, inrange(1) - plus, r, g, b) - call set_pixel(img, out_, inrange(2) - plus, r, g, b) + data(:, out_ + 1, inrange(1) + plus + 1) = rgb + data(:, out_ + 1, inrange(2) + plus + 1) = rgb + data(:, out_ + 1, inrange(1) - plus + 1) = rgb + data(:, out_ + 1, inrange(2) - plus + 1) = rgb end do end do do in_ = inrange(1), inrange(2) do plus = 0, pl % meshlines_width - call set_pixel(img, outrange(1) + plus, in_, r, g, b) - call set_pixel(img, outrange(2) + plus, in_, r, g, b) - call set_pixel(img, outrange(1) - plus, in_, r, g, b) - call set_pixel(img, outrange(2) - plus, in_, r, g, b) + data(:, outrange(1) + plus + 1, in_ + 1) = rgb + data(:, outrange(2) + plus + 1, in_ + 1) = rgb + data(:, outrange(1) - plus + 1, in_ + 1) = rgb + data(:, outrange(2) - plus + 1, in_ + 1) = rgb end do end do @@ -298,12 +296,12 @@ contains ! OUTPUT_PPM writes out a previously generated image to a PPM file !=============================================================================== - subroutine output_ppm(pl, img) + subroutine output_ppm(pl, data) type(ObjectPlot), intent(in) :: pl - type(Image), intent(in) :: img + integer, intent(in) :: data(:,:,:) - integer :: i ! loop index for height - integer :: j ! loop index for width + integer :: y ! loop index for height + integer :: x ! loop index for width integer :: unit_plot ! Open PPM file for writing @@ -311,14 +309,14 @@ contains ! Write header write(unit_plot, '(A2)') 'P6' - write(unit_plot, '(I0,'' '',I0)') img % width, img % height + write(unit_plot, '(I0,'' '',I0)') pl % pixels(1), pl % pixels(2) write(unit_plot, '(A)') '255' ! Write color for each pixel - 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)) + do y = 1, pl % pixels(2) + do x = 1, pl % pixels(1) + write(unit_plot, '(3A1)', advance='no') achar(data(RED, x, y)), & + achar(data(GREEN, x, y)), achar(data(BLUE, x, y)) end do end do @@ -327,7 +325,7 @@ contains end subroutine output_ppm !=============================================================================== -! CREATE_3D_DUMP outputs a binary file that can be input into silomesh for 3D +! CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D ! geometry visualization. It works the same way as create_ppm by dragging a ! particle across the geometry for the specified number of voxels. The first ! 3 int(4)'s in the binary are the number of x, y, and z voxels. The next 3 @@ -338,7 +336,7 @@ contains ! id. For 1 million voxels this produces a file of approximately 15MB. !=============================================================================== - subroutine create_3d_dump(pl) + subroutine create_voxel(pl) type(ObjectPlot), intent(in) :: pl integer :: x, y, z ! voxel location indices @@ -348,7 +346,7 @@ contains integer, target :: data(pl%pixels(3),pl%pixels(2)) integer(HID_T) :: file_id integer(HID_T) :: dspace - integeR(HID_T) :: memspace + integer(HID_T) :: memspace integer(HID_T) :: dset integer(HSIZE_T) :: dims(3) integer(HSIZE_T) :: dims_slab(3) @@ -445,6 +443,6 @@ contains call h5sclose_f(memspace, hdf5_err) call file_close(file_id) - end subroutine create_3d_dump + end subroutine create_voxel end module plot diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 deleted file mode 100644 index ebd8105bef..0000000000 --- a/src/ppmlib.F90 +++ /dev/null @@ -1,127 +0,0 @@ -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 - -!=============================================================================== -! INIT_IMAGE initializes the Image derived type -!=============================================================================== - - 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 - -!=============================================================================== -! ALLOCATE_IMAGE sets the width and height of an image and allocates color -! arrays -!=============================================================================== - - subroutine allocate_image(img, w, h) - - type(Image), intent(inout) :: img - integer, intent(in) :: w ! width of image - integer, intent(in) :: h ! height of image - - ! allocate red, green, and blue array - allocate(img % red(w, h)) - allocate(img % green(w, h)) - allocate(img % blue(w, h)) - - ! set width and height - img % width = w - img % height = h - - end subroutine allocate_image - -!=============================================================================== -! DEALLOCATE_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 - -!=============================================================================== -! INSIDE_IMAGE determines whether a point (x,y) is inside the 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 - -!=============================================================================== -! VALID_IMAGE checks whether the image has a width and height and if its color -! arrays are allocated -!=============================================================================== - - 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 - -!=============================================================================== -! SET_PIXEL sets the colors for a given pixel -!=============================================================================== - - subroutine set_pixel(img, x, y, r, g, b) - - type(Image), intent(inout) :: img - integer, intent(in) :: x, y ! coordinates - integer, intent(in) :: r, g, b ! red, green, and blue - - 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 diff --git a/src/ppmlib.LICENSE b/src/ppmlib.LICENSE deleted file mode 100644 index a733e8d3e9..0000000000 --- a/src/ppmlib.LICENSE +++ /dev/null @@ -1,4 +0,0 @@ -The source code in ppmlib.F90 is adapted from code found on Rosetta Code - which was authored -by Mauro Panigada. The authors of OpenMC have obtained permission from Mauro -Panigada to use and distribution this source code as part of OpenMC.