diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index cf403e8da3..09aea601bb 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1261,7 +1261,7 @@ attributes or sub-elements. These are not used in "voxel" plots: Specifies the RGB color of the regions where no OpenMC cell can be found. Should be three integers separated by spaces. - *Default*: 0 0 0 (white) + *Default*: 0 0 0 (black) :col_spec: Any number of this optional tag may be included in each ```` element, @@ -1301,6 +1301,35 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None + :meshlines: + The ``meshlines`` sub-element allows for plotting the boundaries of + a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per + ``plot`` element, and it must contain as attributes or sub-elements a mesh + type and a linewidth. Optionally, a color may be specified for the overlay: + + :meshtype: + The type of the mesh to be plotted. Valid options are "tally", "entropy", + "ufs", and "cmfd". If plotting "tally" meshes, the id of the mesh to plot + must be specified with the ``id`` sub-element. + + :id: + A single integer id number for the mesh specified on ``tallies.xml`` that + should be plotted. This element is only required for ``meshtype="tally"``. + + :linewidth: + A single integer number of pixels of linewidth to specify for the mesh + boundaries. Specifying this as 0 indicates that lines will be 1 pixel + thick, specifying 1 indicates 3 pixels thick, specifying 2 indicates + 5 pixels thick, etc. + + :color: + Specifies the custom color for the meshlines boundaries. Should be 3 + integers separated by whitespace. This element is optional. + + *Default*: 0 0 0 (black) + + *Default*: None + .. _usersguide_cmfd: ------------------------------ diff --git a/examples/lattice/simple/plots.xml b/examples/lattice/simple/plots.xml index a8c6fbb7c5..a25a5a8280 100644 --- a/examples/lattice/simple/plots.xml +++ b/examples/lattice/simple/plots.xml @@ -5,6 +5,7 @@ 0. 0. 0. 4.0 4.0 400 400 + diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 37edb922c4..cd95568849 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -536,12 +536,12 @@ contains m % n_dimension = 3 allocate(m % dimension(3)) m % dimension = n + + ! determine width + m % width = (m % upper_right - m % lower_left) / m % dimension + end if - ! allocate and determine width - allocate(m % width(3)) - m % width = (m % upper_right - m % lower_left) / m % dimension - ! allocate p allocate(entropy_p(1, m % dimension(1), m % dimension(2), & m % dimension(3))) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 869f9fd63f..e84f072ea8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -36,7 +36,7 @@ contains if (run_mode /= MODE_PLOTTING) call read_cross_sections_xml() call read_geometry_xml() call read_materials_xml() - if (run_mode /= MODE_PLOTTING) call read_tallies_xml() + call read_tallies_xml() if (cmfd_run) call configure_cmfd() end subroutine read_input_xml @@ -517,6 +517,7 @@ contains allocate(entropy_mesh) allocate(entropy_mesh % lower_left(3)) allocate(entropy_mesh % upper_right(3)) + allocate(entropy_mesh % width(3)) ! Copy values call get_node_array(node_entropy, "lower_left", & @@ -548,6 +549,11 @@ contains ! Copy dimensions call get_node_array(node_entropy, "dimension", entropy_mesh % dimension) + + ! Calculate width + entropy_mesh % width = (entropy_mesh % upper_right - & + entropy_mesh % lower_left) / entropy_mesh % dimension + end if ! Turn on Shannon entropy calculation @@ -1895,7 +1901,7 @@ contains end if ! Allocate tally array - if (n_user_tallies > 0) then + if (n_user_tallies > 0 .and. run_mode /= MODE_PLOTTING) then call add_tallies("user", n_user_tallies) end if @@ -2043,6 +2049,9 @@ contains call mesh_dict % add_key(m % id, i) end do + ! We only need the mesh info for plotting + if (run_mode == MODE_PLOTTING) return + ! ========================================================================== ! READ TALLY DATA @@ -2767,20 +2776,25 @@ contains subroutine read_plots_xml() - integer i, j - integer n_cols, col_id, n_comp, n_masks + integer :: i, j + integer :: n_cols, col_id, n_comp, n_masks, n_meshlines + integer :: meshid + integer :: i_mesh integer, allocatable :: iarray(:) logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml character(MAX_LINE_LEN) :: temp_str + character(MAX_WORD_LEN) :: meshtype type(ObjectPlot), pointer :: pl => null() type(Node), pointer :: doc => null() type(Node), pointer :: node_plot => null() type(Node), pointer :: node_col => null() type(Node), pointer :: node_mask => null() + type(Node), pointer :: node_meshlines => null() type(NodeList), pointer :: node_plot_list => null() type(NodeList), pointer :: node_col_list => null() type(NodeList), pointer :: node_mask_list => null() + type(NodeList), pointer :: node_meshline_list => null() ! Check if plots.xml exists filename = trim(path_input) // "plots.xml" @@ -3032,6 +3046,144 @@ contains end do end if + ! Deal with meshlines + call get_node_list(node_plot, "meshlines", node_meshline_list) + n_meshlines = get_list_size(node_meshline_list) + if (n_meshlines /= 0) then + + if (pl % type == PLOT_TYPE_VOXEL) then + message = "Meshlines ignored in voxel plot " // & + trim(to_str(pl % id)) + call warning() + end if + + select case(n_meshlines) + case (0) + ! Skip if no meshlines are specified + case (1) + + ! Get pointer to meshlines + call get_list_item(node_meshline_list, 1, node_meshlines) + + ! Check mesh type + if (check_for_node(node_meshlines, "meshtype")) then + call get_node_value(node_meshlines, "meshtype", meshtype) + else + message = "Must specify a meshtype for meshlines " // & + "specification in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + ! Ensure that there is a linewidth for this meshlines specification + if (check_for_node(node_meshlines, "linewidth")) then + call get_node_value(node_meshlines, "linewidth", & + pl % meshlines_width) + else + message = "Must specify a linewidth for meshlines " // & + "specification in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + ! Check for color + if (check_for_node(node_meshlines, "color")) then + + ! Check and make sure 3 values are specified for RGB + if (get_arraysize_double(node_meshlines, "color") /= 3) then + message = "Bad RGB for meshlines color " // & + "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + call get_node_array(node_meshlines, "color", & + pl % meshlines_color % rgb) + else + + pl % meshlines_color % rgb = (/ 0, 0, 0 /) + + end if + + ! Set mesh based on type + select case (trim(meshtype)) + case ('ufs') + + if (.not. associated(ufs_mesh)) then + message = "No UFS mesh for meshlines on plot " // & + trim(to_str(pl % id)) + call fatal_error() + end if + + pl % meshlines_mesh => ufs_mesh + + case ('cmfd') + + if (.not. cmfd_run) then + message = "Need CMFD run to plot CMFD mesh for meshlines " // & + "on plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + i_mesh = cmfd_tallies(1) % & + filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % & + int_bins(1) + pl % meshlines_mesh => meshes(i_mesh) + + case ('entropy') + + if (.not. associated(entropy_mesh)) then + message = "No entropy mesh for meshlines on plot " // & + trim(to_str(pl % id)) + call fatal_error() + end if + + if (.not. allocated(entropy_mesh % dimension)) then + message = "No dimension specified on entropy mesh for " // & + "meshlines on plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + pl % meshlines_mesh => entropy_mesh + + case ('tally') + + ! Ensure that there is a mesh id if the type is tally + if (check_for_node(node_meshlines, "id")) then + call get_node_value(node_meshlines, "id", meshid) + else + message = "Must specify a mesh id for meshlines tally mesh" // & + "specification in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + ! Check if the specified tally mesh exists + if (mesh_dict % has_key(meshid)) then + pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid)) + if (meshes(meshid) % type /= LATTICE_RECT) then + message = "Non-rectangular mesh specified in meshlines " // & + "for plot " // trim(to_str(pl % id)) + call fatal_error() + end if + else + message = "Could not find mesh " // & + trim(to_str(meshid)) // & + " specified in meshlines for plot " // & + trim(to_str(pl % id)) + call fatal_error() + end if + + case default + message = "Invalid type for meshlines on plot " // & + trim(to_str(pl % id)) // ": " // trim(meshtype) + call fatal_error() + end select + + case default + message = "Mutliple meshlines" // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end select + + end if + ! Deal with masks call get_node_list(node_plot, "mask", node_mask_list) n_masks = get_list_size(node_mask_list) diff --git a/src/plot.F90 b/src/plot.F90 index e7642314f9..91b2520eb4 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -5,6 +5,7 @@ module plot use geometry, only: find_cell, check_cell_overlap use geometry_header, only: Cell, BASE_UNIVERSE use global + use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: deallocate_coord, Particle use plot_header @@ -118,27 +119,24 @@ contains call init_image(img) call allocate_image(img, pl % pixels(1), pl % pixels(2)) + in_pixel = pl % width(1)/dble(pl % pixels(1)) + out_pixel = pl % width(2)/dble(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(2)/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(2) / 2.0 else if (pl % basis == PLOT_BASIS_YZ) then in_i = 2 out_i = 3 - in_pixel = pl % width(1)/dble(pl % pixels(1)) - out_pixel = pl % width(2)/dble(pl % pixels(2)) xyz(1) = pl % origin(1) xyz(2) = pl % origin(2) - pl % width(1) / 2.0 xyz(3) = pl % origin(3) + pl % width(2) / 2.0 @@ -169,6 +167,9 @@ contains p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel end do + ! Draw tally mesh boundaries on the image if requested + if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, img) + ! Write out the ppm to a file call output_ppm(pl,img) @@ -180,6 +181,111 @@ contains end subroutine create_ppm +!=============================================================================== +! DRAW_MESH_LINES draws mesh line boundaries on an image +!=============================================================================== + subroutine draw_mesh_lines(pl, img) + + type(ObjectPlot), pointer, intent(in) :: pl + type(Image), intent(inout) :: img + + logical :: in_mesh + integer :: out_, in_ ! pixel location + integer :: r, g, b ! RGB color for meshlines pixels + integer :: outrange(2), inrange(2) ! range of pixel locations + integer :: i, j ! loop indices + integer :: plus + integer :: ijk_ll(3) ! mesh bin ijk indicies of plot lower left + integer :: ijk_ur(3) ! mesh bin ijk indicies of plot upper right + integer :: outer, inner + real(8) :: frac + real(8) :: width(3) ! real widths of the plot + real(8) :: xyz_ll_plot(3) ! lower left xyz of plot image + real(8) :: xyz_ur_plot(3) ! upper right xyz of plot image + real(8) :: xyz_ll(3) ! lower left xyz + real(8) :: xyz_ur(3) ! upper right xyz + type(StructuredMesh), pointer :: m => null() + + m => pl % meshlines_mesh + + r = pl % meshlines_color % rgb(1) + g = pl % meshlines_color % rgb(2) + b = pl % meshlines_color % rgb(3) + + select case (pl % basis) + case(PLOT_BASIS_XY) + outer = 1 + inner = 2 + case(PLOT_BASIS_XZ) + outer = 1 + inner = 3 + case(PLOT_BASIS_YZ) + outer = 2 + inner = 3 + end select + + xyz_ll_plot = pl % origin + xyz_ur_plot = pl % origin + + xyz_ll_plot(outer) = pl % origin(1) - pl % width(1) / 2.0 + xyz_ll_plot(inner) = pl % origin(2) - pl % width(2) / 2.0 + xyz_ur_plot(outer) = pl % origin(1) + pl % width(1) / 2.0 + xyz_ur_plot(inner) = pl % origin(2) + pl % width(2) / 2.0 + + width = xyz_ur_plot - xyz_ll_plot + + call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) + call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) + + ! sweep through all meshbins on this plane and draw borders + do i = ijk_ll(outer), ijk_ur(outer) + do j = ijk_ll(inner), ijk_ur(inner) + ! check if we're in the mesh for this ijk + if (i > 0 .and. i <= m % dimension(outer) .and. & + j > 0 .and. j <= m % dimension(inner)) then + + ! get xyz's of lower left and upper right of this mesh cell + xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) + xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1) + xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i + xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j + + ! 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)) + frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) + outrange(2) = int(frac * real(img % width, 8)) + + frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(1) = int(frac * real(img % height, 8)) + frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(2) = int(frac * real(img % height, 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) + 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) + end do + end do + + end if + end do + end do + + end subroutine draw_mesh_lines + !=============================================================================== ! OUTPUT_PPM writes out a previously generated image to a PPM file !=============================================================================== diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 6956bfaf93..68eb89a11a 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 @@ -25,6 +26,9 @@ module plot_header real(8) :: width(3) ! xyz widths of plot integer :: basis ! direction of plot slice integer :: pixels(3) ! pixel width/height of plot slice + integer :: meshlines_width ! pixel width of meshlines + type(StructuredMesh), pointer :: meshlines_mesh => null() ! mesh to plot + type(ObjectColor) :: meshlines_color ! Color for meshlines type(ObjectColor) :: not_found ! color for positions where no cell found type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats end type ObjectPlot diff --git a/src/relaxng/plots.rnc b/src/relaxng/plots.rnc index 3c29ba4e0b..a7c657864d 100644 --- a/src/relaxng/plots.rnc +++ b/src/relaxng/plots.rnc @@ -26,6 +26,14 @@ element plots { attribute components { list { xsd:int+ } }) & (element background { list { xsd:int+ } } | attribute background { list { xsd:int+ } }) + }* & + element meshlines { + (element meshtype { ( "tally" | "entropy" | "ufs" | "cmfd" ) } | + attribute meshtype { ( "tally" | "entropy" | "ufs" | "cmfd" ) }) & + (element id { xsd:int } | attribute id { xsd:int })? & + (element linewidth { xsd:int } | attribute linewidth { xsd:int }) & + (element color { list { xsd:int+ } } | + attribute color { list { xsd:int+ } })? }* }* }