From 6a2d8bc0fbb91c816ba8339ae568c431333b7c97 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 7 Mar 2012 10:33:04 -0500 Subject: [PATCH 01/42] fixed plots appearing upside down --- src/plot.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/plot.F90 b/src/plot.F90 index cd2aec0caa..2f82b9dd8f 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -70,7 +70,7 @@ contains 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(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 @@ -79,7 +79,7 @@ contains 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 + xyz(3) = pl % origin(3) + pl % width(2) / 2.0 else if (pl % basis == PLOT_BASIS_YZ) then in_i = 2 out_i = 3 @@ -87,7 +87,7 @@ contains 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 + xyz(3) = pl % origin(3) + pl % width(2) / 2.0 end if ! allocate and initialize particle @@ -132,7 +132,7 @@ contains end do p % coord0 % xyz(in_i) = xyz(in_i) - p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel + p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel end do call output_ppm(pl,img) From 1a9892144cab3eb0a2c887f626658afd00591258 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 7 Mar 2012 11:17:43 -0500 Subject: [PATCH 02/42] added ability to specify colors for specific cells/mats in plots --- src/input_xml.F90 | 55 +++++++++++++++++++++++---- src/xml-fortran/templates/plots_t.xml | 6 +++ 2 files changed, 53 insertions(+), 8 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index fad53be3ea..9b806ade83 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1111,10 +1111,11 @@ contains use xml_data_plots_t - integer i + integer i, j + integer n_cols, col_id logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml - type(Plot), pointer :: pl => null() + type(Plot), pointer :: pl => null() ! Check if plots.xml exists filename = trim(path_input) // "plots.xml" @@ -1145,7 +1146,7 @@ contains if (size(plot_(i) % pixels) == 2) then pl % pixels = plot_(i) % pixels else - message = " must be length 2 in plot " // to_str(i) + message = " must be length 2 in plot " // to_str(pl % id) call fatal_error() end if @@ -1156,7 +1157,7 @@ contains pl % color = PLOT_COLOR_MATS case default message = "Unsupported plot color '" // plot_(i) % color & - // "' in plot " // trim(to_str(i)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -1167,7 +1168,7 @@ contains ! pl % type = PLOT_TYPE_POINTS case default message = "Unsupported plot type '" // plot_(i) % type & - // "' in plot " // trim(to_str(i)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -1180,7 +1181,7 @@ contains pl % basis = PLOT_BASIS_YZ case default message = "Unsupported plot basis '" // plot_(i) % basis & - // "' in plot " // trim(to_str(i)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -1189,7 +1190,7 @@ contains pl % origin = plot_(i) % origin else message = "Origin must be length 3 " & - // "in plot " // trim(to_str(i)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -1201,10 +1202,48 @@ contains pl % width(2) = plot_(i) % width(2) else message = "Bad plot width " & - // "in plot " // trim(to_str(i)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if + ! Copy user specified colors + n_cols = size(plot_(i) % col_spec_) + do j = 1, n_cols + if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then + message = "Bad RGB " & + // "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + col_id = plot_(i) % col_spec_(j) % id + + if (pl % color == PLOT_COLOR_CELLS) then + + if (dict_has_key(cell_dict, col_id)) then + cells(col_id) % rgb(1) = plot_(i) % col_spec_(j) % rgb(1) + cells(col_id) % rgb(2) = plot_(i) % col_spec_(j) % rgb(2) + cells(col_id) % rgb(3) = plot_(i) % col_spec_(j) % rgb(3) + else + message = "Could not find cell " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + else if (pl % color == PLOT_COLOR_MATS) then + + if (dict_has_key(material_dict, col_id)) then + materials(col_id) % rgb(1) = plot_(i) % col_spec_(j) % rgb(1) + materials(col_id) % rgb(2) = plot_(i) % col_spec_(j) % rgb(2) + materials(col_id) % rgb(3) = plot_(i) % col_spec_(j) % rgb(3) + else + message = "Could not find material " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + end if + + end do end do diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml index af8f2c6e27..f372411922 100644 --- a/src/xml-fortran/templates/plots_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -3,6 +3,11 @@ + + + + + @@ -12,6 +17,7 @@ + From c639961c9fa1b1a8527c30dc7085b9744433ffd5 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 15:20:23 -0500 Subject: [PATCH 03/42] Reorganized plotting so colors are stored in Plot objects rather than Cells/Materials, added ability to specify background color, added ability to specify plotting of only a subset of cells --- src/geometry_header.F90 | 1 - src/input_xml.F90 | 93 +++++++++++++++++---------- src/material_header.F90 | 2 - src/output.F90 | 4 +- src/plot.F90 | 24 +++---- src/plot_header.F90 | 14 +++- src/xml-fortran/templates/plots_t.xml | 7 ++ 7 files changed, 93 insertions(+), 52 deletions(-) diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index aa1fcca5f6..fbcf0353a9 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -64,7 +64,6 @@ module geometry_header & surfaces(:) ! List of surfaces bounding cell -- note that ! parentheses, union, etc operators will be listed ! here too - integer, allocatable :: rgb(:) ! plotting color end type Cell ! array index of universe 0 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b806ade83..36ff24fd36 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -261,14 +261,6 @@ contains c % material = cell_(i) % material c % fill = cell_(i) % fill - ! Set plot color - if (plotting) then - allocate(c % rgb(3)) - c % rgb(1) = prn()*255 - c % rgb(2) = prn()*255 - c % rgb(3) = prn()*255 - end if - ! Check to make sure that either material or fill was specified if (c % material == 0 .and. c % fill == 0) then message = "Neither material nor fill was specified for cell " // & @@ -539,14 +531,6 @@ contains ! Copy material id m % id = material_(i) % id - ! Set plot color - if (plotting) then - allocate(m % rgb(3)) - m % rgb(1) = prn()*255 - m % rgb(2) = prn()*255 - m % rgb(3) = prn()*255 - end if - ! Copy density -- the default value for the units is given in the ! material_t.xml file and doesn't need to be specified here, hence case ! default results in an error. @@ -1143,6 +1127,7 @@ contains pl % id = plot_(i) % id pl % aspect = plot_(i) % aspect + ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then pl % pixels = plot_(i) % pixels else @@ -1150,17 +1135,18 @@ contains call fatal_error() end if - select case (plot_(i) % color) - case ("cell") - pl % color = PLOT_COLOR_CELLS - case ("mat") - pl % color = PLOT_COLOR_MATS - case default - message = "Unsupported plot color '" // plot_(i) % color & - // "' in plot " // trim(to_str(pl % id)) + ! Copy plot background color + if (size(plot_(i) % background) == 3) then + pl % not_found % rgb = plot_(i) % background + else if (size(plot_(i) % background) == 0) then + pl % not_found % rgb = (/ 0, 0, 0 /) + else + message = "Bad background RGB " & + // "in plot " // trim(to_str(pl % id)) call fatal_error() - end select + end if + ! Copy plot type select case (plot_(i) % type) case ("slice") pl % type = PLOT_TYPE_SLICE @@ -1172,6 +1158,7 @@ contains call fatal_error() end select + ! Copy plot basis select case (plot_(i) % basis) case ("xy") pl % basis = PLOT_BASIS_XY @@ -1206,6 +1193,34 @@ contains call fatal_error() end if + ! Copy plot color type and initialize all colors randomly + select case (plot_(i) % color) + case ("cell") + + pl % color_by = PLOT_COLOR_CELLS + allocate(pl % colors(n_cells)) + do j = 1, n_cells + pl % colors(j) % rgb(1) = prn()*255 + pl % colors(j) % rgb(2) = prn()*255 + pl % colors(j) % rgb(3) = prn()*255 + end do + + case ("mat") + + pl % color_by = PLOT_COLOR_MATS + allocate(pl % colors(n_materials)) + do j = 1, n_materials + pl % colors(j) % rgb(1) = prn()*255 + pl % colors(j) % rgb(2) = prn()*255 + pl % colors(j) % rgb(3) = prn()*255 + end do + + case default + message = "Unsupported plot color type '" // plot_(i) % color & + // "' in plot " // trim(to_str(pl % id)) + call fatal_error() + end select + ! Copy user specified colors n_cols = size(plot_(i) % col_spec_) do j = 1, n_cols @@ -1217,24 +1232,20 @@ contains col_id = plot_(i) % col_spec_(j) % id - if (pl % color == PLOT_COLOR_CELLS) then + if (pl % color_by == PLOT_COLOR_CELLS) then if (dict_has_key(cell_dict, col_id)) then - cells(col_id) % rgb(1) = plot_(i) % col_spec_(j) % rgb(1) - cells(col_id) % rgb(2) = plot_(i) % col_spec_(j) % rgb(2) - cells(col_id) % rgb(3) = plot_(i) % col_spec_(j) % rgb(3) + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find cell " // trim(to_str(col_id)) // & " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if - else if (pl % color == PLOT_COLOR_MATS) then + else if (pl % color_by == PLOT_COLOR_MATS) then if (dict_has_key(material_dict, col_id)) then - materials(col_id) % rgb(1) = plot_(i) % col_spec_(j) % rgb(1) - materials(col_id) % rgb(2) = plot_(i) % col_spec_(j) % rgb(2) - materials(col_id) % rgb(3) = plot_(i) % col_spec_(j) % rgb(3) + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find material " // trim(to_str(col_id)) // & " specified in plot " // trim(to_str(pl % id)) @@ -1245,6 +1256,22 @@ contains end do + ! Alter colors based on mask information + if (size(plot_(i) % mask_) /= 1) then + message = "Mutliple masks" // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + else + do j=1,size(pl % colors) + if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then + pl % colors(j) % rgb = plot_(i) % mask_(1) % background + end if + end do + end if + + + + end do end subroutine read_plots_xml diff --git a/src/material_header.F90 b/src/material_header.F90 index d8e9ca4ea2..916578923b 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -21,8 +21,6 @@ module material_header character(12) :: sab_name ! name of S(a,b) table integer :: sab_table = 0 ! index in sab_tables integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table - - integer, allocatable :: rgb(:) ! plotting color end type Material end module material_header diff --git a/src/output.F90 b/src/output.F90 index c4209ace13..0876d266c2 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -914,7 +914,7 @@ contains write(ou,100) "Width:", trim(to_str(pl % width(1))) // & " " // trim(to_str(pl % width(2))) - write(ou,100) "Coloring:", trim(to_str(pl % color)) + write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & trim(to_str(pl % pixels(2))) @@ -924,7 +924,7 @@ contains write(ou,100) "Width:", trim(to_str(pl % width(1))) // & " " // trim(to_str(pl % width(2))) // " " & // trim(to_str(pl % width(3))) - write(ou,100) "Coloring:", trim(to_str(pl % color)) + write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect)) end if diff --git a/src/plot.F90 b/src/plot.F90 index 2f82b9dd8f..ef93fecf4f 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -106,19 +106,19 @@ contains call find_cell(found_cell) if (.not. found_cell) then - r = 255 - g = 255 - b = 255 + r = pl % not_found % rgb(1) + g = pl % not_found % rgb(2) + b = pl % not_found % rgb(3) else - c => cells(p % coord % cell) - if (pl % color == PLOT_COLOR_MATS) then - r = materials(c % material) % rgb(1) - g = materials(c % material) % rgb(2) - b = materials(c % material) % rgb(3) - else if (pl % color == PLOT_COLOR_CELLS) then - r = c % rgb(1) - g = c % rgb(2) - b = c % rgb(3) + if (pl % color_by == PLOT_COLOR_MATS) then + c => cells(p % coord % cell) + r = pl % colors(c % material) % rgb(1) + g = pl % colors(c % material) % rgb(2) + b = pl % colors(c % material) % rgb(3) + else if (pl % color_by == PLOT_COLOR_CELLS) then + r = pl % colors(p % coord % cell) % rgb(1) + g = pl % colors(p % coord % cell) % rgb(2) + b = pl % colors(p % coord % cell) % rgb(3) else r = 0 g = 0 diff --git a/src/plot_header.F90 b/src/plot_header.F90 index bfe6f9b071..b2f849ea99 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -3,18 +3,28 @@ module plot_header implicit none !=============================================================================== -! PLOT hold plot information +! ObjectColor holds color information for plotted objects +!=============================================================================== + + type ObjectColor + integer :: rgb(3) + end type ObjectColor + +!=============================================================================== +! PLOT holds plot information !=============================================================================== type Plot integer :: id ! Unique ID integer :: type ! Type - integer :: color ! quantity to color regions by + integer :: color_by ! quantity to color regions by real(8) :: origin(3) ! xyz center of plot location real(8) :: aspect ! spacing between rays in raytracer real(8) :: width(3) ! xyz widths of plot integer :: basis ! direction of plot slice integer :: pixels(2) ! 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 end type Plot integer :: PLOT_TYPE_SLICE = 1 diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml index f372411922..e3f6903fe6 100644 --- a/src/xml-fortran/templates/plots_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -8,6 +8,11 @@ + + + + + @@ -17,7 +22,9 @@ + + From 41b40e58dd3cd28ea0c6d55c202ca05ed0dc8171 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 15:49:12 -0500 Subject: [PATCH 04/42] added ability to specify filename for plot output ppm --- src/DEPENDENCIES | 2 ++ src/input_xml.F90 | 4 ++++ src/plot.F90 | 4 +--- src/plot_header.F90 | 3 +++ src/xml-fortran/templates/plots_t.xml | 1 + 5 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 72f98c217d..f2cb4f5e12 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -180,6 +180,8 @@ physics.o: search.o physics.o: string.o physics.o: tally.o +plot_header.o: constants.o + plot.o: constants.o plot.o: error.o plot.o: geometry.o diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 36ff24fd36..3a7bb02806 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1127,6 +1127,10 @@ contains pl % id = plot_(i) % id pl % aspect = plot_(i) % aspect + ! Set output file path + pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & + "_" // trim(plot_(i) % filename) // ".ppm" + ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then pl % pixels = plot_(i) % pixels diff --git a/src/plot.F90 b/src/plot.F90 index ef93fecf4f..61d0ea74dd 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -151,10 +151,8 @@ contains type(Image), intent(in) :: img integer :: i, j - character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file - path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm" - open(UNIT=UNIT_PLOT, FILE=path_plot) + open(UNIT=UNIT_PLOT, FILE=pl % path_plot) write(UNIT_PLOT, '(A2)') 'P6' write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height diff --git a/src/plot_header.F90 b/src/plot_header.F90 index b2f849ea99..294a7d0f5d 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -1,5 +1,7 @@ module plot_header + use constants + implicit none !=============================================================================== @@ -16,6 +18,7 @@ module plot_header type Plot integer :: id ! Unique ID + character(MAX_LINE_LEN) :: path_plot ! path for plot file integer :: type ! Type integer :: color_by ! quantity to color regions by real(8) :: origin(3) ! xyz center of plot location diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml index e3f6903fe6..f01a21a836 100644 --- a/src/xml-fortran/templates/plots_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -15,6 +15,7 @@ + From 0dd4b878458e1d1a6e3b655da34a59832b45056c Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 16:12:35 -0500 Subject: [PATCH 05/42] removed unused plot item 'aspect' --- src/input_xml.F90 | 1 - src/output.F90 | 9 --------- src/plot_header.F90 | 1 - src/xml-fortran/templates/plots_t.xml | 1 - 4 files changed, 12 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3a7bb02806..b19c5b4a86 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1125,7 +1125,6 @@ contains ! Copy data into plots pl % id = plot_(i) % id - pl % aspect = plot_(i) % aspect ! Set output file path pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & diff --git a/src/output.F90 b/src/output.F90 index 0876d266c2..1a9141c0e0 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -918,15 +918,6 @@ contains write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & trim(to_str(pl % pixels(2))) - - else if (pl % type == PLOT_TYPE_POINTS) then - - write(ou,100) "Width:", trim(to_str(pl % width(1))) // & - " " // trim(to_str(pl % width(2))) // " " & - // trim(to_str(pl % width(3))) - write(ou,100) "Coloring:", trim(to_str(pl % color_by)) - write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect)) - end if write(ou,*) diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 294a7d0f5d..c0ae470643 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -22,7 +22,6 @@ module plot_header integer :: type ! Type integer :: color_by ! quantity to color regions by real(8) :: origin(3) ! xyz center of plot location - real(8) :: aspect ! spacing between rays in raytracer real(8) :: width(3) ! xyz widths of plot integer :: basis ! direction of plot slice integer :: pixels(2) ! pixel width/height of plot slice diff --git a/src/xml-fortran/templates/plots_t.xml b/src/xml-fortran/templates/plots_t.xml index f01a21a836..b3b86fc44e 100644 --- a/src/xml-fortran/templates/plots_t.xml +++ b/src/xml-fortran/templates/plots_t.xml @@ -22,7 +22,6 @@ - From be27d56c9421a499a6ece8b76adeecaf662af6ea Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 16:39:23 -0500 Subject: [PATCH 06/42] updated plotting documentation to describe color specification and masking --- docs/source/usersguide/input.rst | 123 +++++++++++++++++++++++++------ 1 file changed, 99 insertions(+), 24 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 89431ba114..e213d65c48 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -573,40 +573,115 @@ tallies. This element should be followed by "yes" or "no" Geometry Plotting Specification -- plot.xml ------------------------------------------- -A rudimentary plotting capability is available in OpenMC by specifying a -plot.xml file and subsequently running with the command-line flag ``-plot``. The -root element of the plot.xml is simply ```` and four sub-elements can be -defined to configure the plotting range and resolution. +A basic 2D plotting capability is available in OpenMC by creating a +plots.xml file and subsequently running with the command-line flag ``-plot``. The +root element of the plot.xml is simply ```` and any number output +figures can be defined with ```` sub-elements. -``origin`` Element +``plot`` Element ------------------ -The ``origin`` element has no attributes/sub-elements and indicates the -Cartesian coordinates of the center of the plot. +Each plot must contain a combination of the following attributes or sub-elements: - *Default*: None + :id: + The unique ``id`` of the plot. -``width`` Element ------------------ + *Default*: None - Required entry -The ``width`` element has no attributes/sub-elements and indicates the width of -the plot in each of the basis directions. + :filename: + Filename for the output plot file. - *Default*: None + *Default*: "plot" -``basis`` Element ------------------ + :color: + Keyword for plot coloring. This can only be either ``cell`` or ``mat``, + which colors regions by cells and materials, respectively. -The ``basis`` element has no attributes/sub-elements and indicates the specified -basis for plotting. The only option option currently accepted is "xy". + *Default*: ``cell`` - *Default*: xy + :origin: + Specifies the XYZ coordinate of the center of the plot. Should be 3 floats + separated by spaces. -``pixel`` Element ------------------ + *Default*: None - Required entry -The ``pixel`` element has no attributes/sub-elements and indicates the distance -between horizontal rays sent through the geometry to record surface crossings. A -smaller ``pixel`` will result in a higher-resolution plot. + :width: + Specifies the width of the plot along each of the basis directions. + Should be 2 or 3 floats separated by spaces for 2D plots and 3D plots, + respectively. - *Default*: 0.01 + *Default*: None - Required entry + + :type: + Keyword for type of plot to be produced. Currently only ``slice`` plots are + implemented, which create 2D pixel maps saved in the PPM file format. PPM + files can be displayed in most viewers (e.g. the default Gnome viewer, + IrfanView, etc.). + + .. note:: Since the PPM format is saved without any kind of compression, + the resulting file sizes can be quite large. Saving the image in + the PNG format can often times reduce the file size by orders of + magnitude without any loss of image quality. + + *Default*: "slice" + +``plot`` elements of ``type`` ``slice`` also contain the following attributes or +sub-elements: + + :basis: + Keyword specifying the plane of the plot for ``slice`` type plots. Can be + one of: ``xy``, ``xz``, ``yz``. + + *Default*: ``xy`` + + :pixels: + Specifies the number of pixes to be used along each of the basis directions for + ``slice`` plots. Should be 2 integers separated by spaces. + + .. warning:: The ``pixels`` input determines the output file size. For the PPM + format, 10 million pixels will result in a file just under 30 MB in + size. + + .. warning:: If the aspect ratio defined in ``pixels`` does not match the aspect + ratio defined in ``width`` the plot may appear stretched or squeezed. + + .. warning:: Geometry features along a basis direction smaller than ``width``/``pixels`` + along that basis direction may not appear in the plot. + + *Default*: None - Required entry for ``slice`` plots + + :background: + Specifies the RGB color of the regions where no OpenMC cell can be found. Should + be 3 integers deparated by spaces. + + *Default*: 0 0 0 (white) + + :col_spec: + Any number of this optional tag may be included in each ``plot`` element, which can + override the default random colors for cells or materials. Each ``col_spec`` + element must contain ``id`` and ``rgb`` sub-elements. + + :id: + Specifies the cell or material unique id for the color specification. + + :rgb: + Specifies the custom color for the cell or material. Should be 3 intergers separated + by spaces. + + *Default*: None + + :mask: + The special ``mask`` sub-element allows for the selective plotting of *only* + user-specified cells or materials. Only one ``mask`` element is allowed per ``plot`` + element, and it must contain as atributes or sub-elements a background masking color and + a list of cells or materials to plot: + + :components: + List of unique ``id`` numbers of the cells or materials to plot. Should be any number + of integers separated by spaces. + + :background: + Color to apply to all cells or materials not in the ``components`` list of cells or + materials to plot. This overrides any ``col_spec`` color specifications. + + *Default*: None From 79b0aae5ffccbc71a2654fda7119f2c3d007c5d8 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 16:55:50 -0500 Subject: [PATCH 07/42] minor bugfix to plot slice xml parsing for mask elements --- src/input_xml.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b19c5b4a86..285e9974e4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1260,11 +1260,11 @@ contains end do ! Alter colors based on mask information - if (size(plot_(i) % mask_) /= 1) then + if (size(plot_(i) % mask_) > 1) then message = "Mutliple masks" // & " specified in plot " // trim(to_str(pl % id)) call fatal_error() - else + else if (.not. size(plot_(i) % mask_) == 0) then do j=1,size(pl % colors) if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then pl % colors(j) % rgb = plot_(i) % mask_(1) % background From 6b638f2d7b4e4631cfba37c98c2e5f86d24f18e0 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 9 Mar 2012 17:05:05 -0500 Subject: [PATCH 08/42] changed default plot background color back to white --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 285e9974e4..3f74d514a6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1142,7 +1142,7 @@ contains if (size(plot_(i) % background) == 3) then pl % not_found % rgb = plot_(i) % background else if (size(plot_(i) % background) == 0) then - pl % not_found % rgb = (/ 0, 0, 0 /) + pl % not_found % rgb = (/ 255, 255, 255 /) else message = "Bad background RGB " & // "in plot " // trim(to_str(pl % id)) From a1b056012e51968afe83418cab2d7f990807528a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Mar 2012 11:13:54 -0400 Subject: [PATCH 09/42] Fixed inputs in examples directory. --- examples/basic/materials.xml | 11 ++++++----- examples/basic/settings.xml | 5 ----- examples/basic/tallies.xml | 6 +++--- examples/lattice/nested/geometry.xml | 4 ++-- examples/lattice/nested/plot.xml | 6 ------ examples/lattice/nested/plots.xml | 10 ++++++++++ examples/lattice/nested/tallies.xml | 6 +++--- examples/lattice/simple/geometry.xml | 2 +- examples/lattice/simple/plot.xml | 6 ------ examples/lattice/simple/plots.xml | 10 ++++++++++ examples/lattice/simple/tallies.xml | 6 +++--- examples/reflective/materials.xml | 5 +++-- examples/reflective/settings.xml | 5 ----- 13 files changed, 41 insertions(+), 41 deletions(-) delete mode 100644 examples/lattice/nested/plot.xml create mode 100644 examples/lattice/nested/plots.xml delete mode 100644 examples/lattice/simple/plot.xml create mode 100644 examples/lattice/simple/plots.xml diff --git a/examples/basic/materials.xml b/examples/basic/materials.xml index 88e20c69ef..e50a423112 100644 --- a/examples/basic/materials.xml +++ b/examples/basic/materials.xml @@ -1,17 +1,18 @@ - + 70c + - + - - - + + + diff --git a/examples/basic/settings.xml b/examples/basic/settings.xml index a9a83dfe54..e9c2e78429 100644 --- a/examples/basic/settings.xml +++ b/examples/basic/settings.xml @@ -1,11 +1,6 @@ - - - /home/paulromano/openmc/cross_sections_serpent.xml - - 15 diff --git a/examples/basic/tallies.xml b/examples/basic/tallies.xml index acafe76742..a434deae11 100644 --- a/examples/basic/tallies.xml +++ b/examples/basic/tallies.xml @@ -5,7 +5,7 @@ 100 - total scatter nu-scatter absorption fission nu-fission + total scatter nu-scatter absorption fission nu-fission @@ -13,7 +13,7 @@ 100 0 20.0 - total scatter nu-scatter absorption fission nu-fission + total scatter nu-scatter absorption fission nu-fission @@ -22,7 +22,7 @@ 0 20.0 0 20.0 - scatter nu-scatter nu-fission + scatter nu-scatter nu-fission diff --git a/examples/lattice/nested/geometry.xml b/examples/lattice/nested/geometry.xml index 4ce62280a2..f223cf3bb4 100644 --- a/examples/lattice/nested/geometry.xml +++ b/examples/lattice/nested/geometry.xml @@ -14,7 +14,7 @@ rectangular 2 2 - -1.0 -1.0 + -1.0 -1.0 1.0 1.0 1 2 @@ -26,7 +26,7 @@ rectangular 2 2 - -2.0 -2.0 + -2.0 -2.0 2.0 2.0 5 5 diff --git a/examples/lattice/nested/plot.xml b/examples/lattice/nested/plot.xml deleted file mode 100644 index d6a0ca3cdf..0000000000 --- a/examples/lattice/nested/plot.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - 0. 0. 0. - 3.999 3.999 - 0.01 - diff --git a/examples/lattice/nested/plots.xml b/examples/lattice/nested/plots.xml new file mode 100644 index 0000000000..a8c6fbb7c5 --- /dev/null +++ b/examples/lattice/nested/plots.xml @@ -0,0 +1,10 @@ + + + + + 0. 0. 0. + 4.0 4.0 + 400 400 + + + diff --git a/examples/lattice/nested/tallies.xml b/examples/lattice/nested/tallies.xml index a389dc905c..271e81fe79 100644 --- a/examples/lattice/nested/tallies.xml +++ b/examples/lattice/nested/tallies.xml @@ -4,7 +4,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 @@ -13,7 +13,7 @@ 1 - total + total - \ No newline at end of file + diff --git a/examples/lattice/simple/geometry.xml b/examples/lattice/simple/geometry.xml index 4a77229392..b5bad90b9c 100644 --- a/examples/lattice/simple/geometry.xml +++ b/examples/lattice/simple/geometry.xml @@ -12,7 +12,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 1 2 1 2 diff --git a/examples/lattice/simple/plot.xml b/examples/lattice/simple/plot.xml deleted file mode 100644 index d6a0ca3cdf..0000000000 --- a/examples/lattice/simple/plot.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - 0. 0. 0. - 3.999 3.999 - 0.01 - diff --git a/examples/lattice/simple/plots.xml b/examples/lattice/simple/plots.xml new file mode 100644 index 0000000000..a8c6fbb7c5 --- /dev/null +++ b/examples/lattice/simple/plots.xml @@ -0,0 +1,10 @@ + + + + + 0. 0. 0. + 4.0 4.0 + 400 400 + + + diff --git a/examples/lattice/simple/tallies.xml b/examples/lattice/simple/tallies.xml index a389dc905c..271e81fe79 100644 --- a/examples/lattice/simple/tallies.xml +++ b/examples/lattice/simple/tallies.xml @@ -4,7 +4,7 @@ rectangular 4 4 - -2.0 -2.0 + -2.0 -2.0 1.0 1.0 @@ -13,7 +13,7 @@ 1 - total + total - \ No newline at end of file + diff --git a/examples/reflective/materials.xml b/examples/reflective/materials.xml index eb55455f7a..3ae574ff3e 100644 --- a/examples/reflective/materials.xml +++ b/examples/reflective/materials.xml @@ -1,10 +1,11 @@ - + 70c + - + diff --git a/examples/reflective/settings.xml b/examples/reflective/settings.xml index 007b9cdd95..02ad736305 100644 --- a/examples/reflective/settings.xml +++ b/examples/reflective/settings.xml @@ -1,11 +1,6 @@ - - - /home/paulromano/openmc/cross_sections_serpent.xml - - 500 From 2662a1b43b77ed266d9c82ad5e52ecdaa10fbe8c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Mar 2012 12:35:26 -0400 Subject: [PATCH 10/42] Added MIT license. --- LICENSE | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/LICENSE b/LICENSE index e69de29bb2..8c94262bd2 100644 --- a/LICENSE +++ b/LICENSE @@ -0,0 +1,18 @@ +Copyright (c) 2011-2012 Massachusetts Institute of Technology + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. From 1e70be6f9b04fd4c2738cc3b8b65af208ddf1c0f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 15 Mar 2012 17:50:01 -0400 Subject: [PATCH 11/42] Added capability to due batching and multiple generations per batch. Not tested yet. --- src/error.F90 | 7 +- src/global.F90 | 17 +++-- src/hdf5_interface.F90 | 17 +++-- src/input_xml.F90 | 17 +++-- src/intercycle.F90 | 84 ++++++++++++------------ src/main.F90 | 77 ++++++++++++---------- src/output.F90 | 13 ++-- src/source.F90 | 5 +- src/tally.F90 | 19 +++--- src/xml-fortran/templates/settings_t.xml | 3 +- 10 files changed, 143 insertions(+), 116 deletions(-) diff --git a/src/error.F90 b/src/error.F90 index 4f4c53d78a..59cc6d80e0 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -76,9 +76,10 @@ contains write(ERROR_UNIT,*) ! Write information on current cycle and particle - if (current_cycle > 0) then - write(ERROR_UNIT,'(1X,A,I11) ') 'Cycle: ', current_cycle - write(ERROR_UNIT,'(1X,A,I11)') 'Particle: ', p % id + if (current_batch > 0) then + write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch + write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen + write(ERROR_UNIT,'(1X,A,I12)') 'Particle: ', p % id write(ERROR_UNIT,*) end if diff --git a/src/global.F90 b/src/global.F90 index bd4b7471c4..d7edf18215 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -125,11 +125,13 @@ module global ! ============================================================================ ! CRITICALITY SIMULATION VARIABLES - integer(8) :: n_particles = 10000 ! # of particles per cycle - integer :: n_cycles = 500 ! # of cycles - integer :: n_inactive = 50 ! # of inactive cycles - integer :: n_active ! # of active cycles - integer :: current_cycle = 0 ! current cycle + integer(8) :: n_particles ! # of particles per generation + integer :: n_batches ! # of batches + integer :: n_inactive ! # of inactive batches + integer :: n_active ! # of active batches + integer :: gen_per_batch = 1 ! # of generations per batch + integer :: current_batch = 0 ! current batch + integer :: current_gen = 0 ! current generation ! External source type(ExtSource), target :: external_source @@ -143,7 +145,7 @@ module global integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - ! cycle keff + ! single-genreation keff real(8) :: keff = ONE real(8) :: keff_std @@ -217,7 +219,8 @@ module global ! Trace for single particle logical :: trace - integer :: trace_cycle + integer :: trace_batch + integer :: trace_gen integer(8) :: trace_particle contains diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index fafe30109c..cfc7383291 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -77,19 +77,22 @@ contains call hdf5_make_double(hdf5_output_file, "n_particles", real(n_particles,8)) ! Use H5LT interface to write n_cycles, n_inactive, and n_active - call hdf5_make_integer(hdf5_output_file, "n_cycles", n_cycles) + call hdf5_make_integer(hdf5_output_file, "n_batches", n_batches) call hdf5_make_integer(hdf5_output_file, "n_inactive", n_inactive) call hdf5_make_integer(hdf5_output_file, "n_active", n_active) + call hdf5_make_integer(hdf5_output_file, "gen_per_batch", gen_per_batch) ! Add description of each variable call h5ltset_attribute_string_f(hdf5_output_file, "n_particles", & "description", "Number of particles per cycle", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_cycles", & - "description", "Total number of cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "n_batches", & + "description", "Total number of batches", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_inactive", & "description", "Number of inactive cycles", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_active", & "description", "Number of active cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "gen_per_batch", & + "description", "Number of generations per batch", hdf5_err) end if call hdf5_write_geometry() @@ -788,7 +791,7 @@ contains call h5ltset_attribute_string_f(timing_group, "time_transport", & "description", "Time in transport only (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_intercycle", & - "description", "Total time between cycles (s)", hdf5_err) + "description", "Total time between generations (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_tallies", & "description", "Time between cycles accumulating tallies (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_sample", & @@ -796,16 +799,16 @@ contains call h5ltset_attribute_string_f(timing_group, "time_sendrecv", & "description", "Time between cycles SEND/RECVing source sites (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_inactive", & - "description", "Total time in inactive cycles (s)", hdf5_err) + "description", "Total time in inactive batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_active", & - "description", "Total time in active cycles (s)", hdf5_err) + "description", "Total time in active batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_finalize", & "description", "Total time for finalization (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_total", & "description", "Total time elapsed (s)", hdf5_err) ! Write calculation rate - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) call hdf5_make_double(timing_group, "neutrons_per_second", speed) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3f74d514a6..29c72309a9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -91,7 +91,7 @@ contains end if ! Criticality information - if (criticality % cycles > 0) then + if (criticality % batches > 0) then problem_type = PROB_CRITICALITY ! Check number of particles @@ -102,9 +102,13 @@ contains n_particles = str_to_int(criticality % particles) ! Copy cycle information - n_cycles = criticality % cycles - n_inactive = criticality % inactive - n_active = n_cycles - n_inactive + n_batches = criticality % batches + n_inactive = criticality % inactive + n_active = n_batches - n_inactive + gen_per_batch = criticality % generations_per_batch + else + message = "Need to specify number of batches with tag." + call fatal_error() end if ! Verbosity @@ -151,8 +155,9 @@ contains ! Particle trace if (associated(trace_)) then - trace_cycle = trace_(1) - trace_particle = trace_(2) + trace_batch = trace_(1) + trace_gen = trace_(2) + trace_particle = trace_(3) end if ! Shannon Entropy mesh diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 5ecd79fe84..6d3ae131e8 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -10,7 +10,7 @@ module intercycle use random_lcg, only: prn, set_particle_seed, prn_skip use search, only: binary_search use string, only: to_str - use tally, only: accumulate_cycle_estimate + use tally, only: accumulate_batch_estimate use tally_header, only: TallyObject use timing, only: timer_start, timer_stop @@ -22,8 +22,8 @@ contains !=============================================================================== ! SYNCHRONIZE_BANK samples source sites from the fission sites that were -! accumulated during the cycle. This routine is what allows this Monte Carlo to -! scale to large numbers of processors where other codes cannot. +! accumulated during the generation. This routine is what allows this Monte +! Carlo to scale to large numbers of processors where other codes cannot. !=============================================================================== subroutine synchronize_bank() @@ -73,7 +73,7 @@ contains total = n_bank #endif - ! If there are not that many particles per cycle, it's possible that no + ! If there are not that many particles per generation, it's possible that no ! fission sites were created at all on a single processor. Rather than add ! extra logic to treat this circumstance, we really want to ensure the user ! runs enough particles to avoid this in the first place. @@ -83,12 +83,12 @@ contains call fatal_error() end if - ! Make sure all processors start at the same point for random sampling by - ! using the current_cycle as the starting seed. Then skip ahead in the - ! sequence using the starting index in the 'global' fission bank for each - ! processor Skip ahead however many random numbers are needed + ! Make sure all processors start at the same point for random sampling. Then + ! skip ahead in the sequence using the starting index in the 'global' + ! fission bank for each processor. - call set_particle_seed(int(current_cycle,8)) + call set_particle_seed(int((current_batch - 1)*gen_per_batch + & + current_gen,8)) call prn_skip(start) ! Determine how many fission sites we need to sample from the source bank @@ -257,20 +257,22 @@ contains ! Send we initiated a series of asynchronous ISENDs and IRECVs, now we have ! to ensure that the data has actually been communicated before moving on to - ! the next cycle + ! the next generation call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err) - ! Deallocate space for bank_position on the last cycle - if (current_cycle == n_cycles) deallocate(bank_position) + ! Deallocate space for bank_position on the very last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(bank_position) #else source_bank = temp_sites(1:n_particles) #endif call timer_stop(time_ic_sendrecv) - ! Deallocate space for the temporary source bank on the last cycle - if (current_cycle == n_cycles) deallocate(temp_sites) + ! Deallocate space for the temporary source bank on the last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(temp_sites) end subroutine synchronize_bank @@ -379,27 +381,22 @@ contains end subroutine shannon_entropy !=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them +! CALCULATE_KEFF calculates the single batch estimate of keff as well as the +! mean and standard deviation of the mean for active batches and displays them !=============================================================================== subroutine calculate_keff() - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff + integer :: n ! active batch number + real(8) :: k_batch ! single batch estimate of keff #ifdef MPI real(8) :: global_temp(N_GLOBAL_TALLIES) #endif - message = "Calculate cycle keff..." + message = "Calculate batch keff..." call write_message(8) - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value #ifdef MPI - global_tallies(K_ANALOG) % value = n_bank * keff - ! Copy global tallies into array to be reduced global_temp = global_tallies(:) % value @@ -416,49 +413,52 @@ contains ! Reset value on other processors global_tallies(:) % value = ZERO end if -#else - global_tallies(K_ANALOG) % value = n_bank * keff #endif ! Collect statistics and print output if (master) then - k_cycle = global_tallies(K_ANALOG) % value/n_particles + k_batch = global_tallies(K_ANALOG) % value/(n_particles*gen_per_batch) - if (current_cycle > n_inactive) then - ! Active cycle number - n = current_cycle - n_inactive + if (tallies_on) then + ! Active batch number + n = current_batch - n_inactive - ! Accumulate single cycle realizations of k - call accumulate_cycle_estimate(global_tallies) + ! Accumulate single batch realizations of k + call accumulate_batch_estimate(global_tallies) ! Determine mean and standard deviation of mean keff = global_tallies(K_ANALOG) % sum/n keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n) - ! Display output for this cycle - if (current_cycle > n_inactive + 1) then + ! Display output for this batch + if (current_batch > n_inactive + 1) then if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=103) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std else - write(UNIT=OUTPUT_UNIT, FMT=101) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=101) current_batch, k_batch, & keff, keff_std end if else if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if end if else - ! Display output for inactive cycle + ! Display output for inactive batch if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if - keff = k_cycle + + ! Set keff + keff = k_batch + + ! Reset tally values + global_tallies(:) % value = ZERO end if end if diff --git a/src/main.F90 b/src/main.F90 index af9824a6a0..31092f79f6 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -11,7 +11,7 @@ program main use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str - use tally, only: synchronize_tallies + use tally, only: synchronize_tallies, add_to_score use timing, only: timer_start, timer_stop #ifdef MPI @@ -54,51 +54,68 @@ contains ! Display column titles if (entropy_on) then - message = " Cycle k(cycle) Entropy Average k" + message = " Batch k(batch) Entropy Average k" call write_message(1) message = " ===== ======== ======= ===================" call write_message(1) else - message = " Cycle k(cycle) Average k" + message = " Batch k(batch) Average k" call write_message(1) message = " ===== ======== ===================" call write_message(1) end if ! ========================================================================== - ! LOOP OVER CYCLES - CYCLE_LOOP: do current_cycle = 1, n_cycles + ! LOOP OVER BATCHES + BATCH_LOOP: do current_batch = 1, n_batches - message = "Simulating cycle " // trim(to_str(current_cycle)) // "..." + message = "Simulating batch " // trim(to_str(current_batch)) // "..." call write_message(8) - ! Set all tallies to zero - n_bank = 0 - ! ======================================================================= - ! LOOP OVER HISTORIES + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch - ! Start timer for transport - call timer_start(time_transport) + ! Set all tallies to zero + n_bank = 0 - HISTORY_LOOP: do i = 1, work + ! ==================================================================== + ! LOOP OVER HISTORIES - ! grab source particle from bank - call get_source_particle(i) + ! Start timer for transport + call timer_start(time_transport) - ! transport particle - call transport() + HISTORY_LOOP: do i = 1, work - end do HISTORY_LOOP + ! grab source particle from bank + call get_source_particle(i) - ! Accumulate time for transport - call timer_stop(time_transport) + ! transport particle + call transport() - ! ======================================================================= - ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + end do HISTORY_LOOP - ! Start timer for inter-cycle synchronization - call timer_start(time_intercycle) + ! Accumulate time for transport + call timer_stop(time_transport) + + ! ==================================================================== + ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + + ! Start timer for inter-cycle synchronization + call timer_start(time_intercycle) + + ! Distribute fission bank across processors evenly + call synchronize_bank() + + ! Add to analog estimate of keff -- since the creation of bank sites + ! was originally weighted by the last cycle keff, we need to multiply + ! by that keff to get the current cycle's value + call add_to_score(global_tallies(K_ANALOG), n_bank * keff) + + ! Stop timer for inter-cycle synchronization + call timer_stop(time_intercycle) + + end do GENERATION_LOOP ! Collect tallies if (tallies_on) then @@ -110,25 +127,17 @@ contains ! Calculate shannon entropy if (entropy_on) call shannon_entropy() - ! Distribute fission bank across processors evenly - call synchronize_bank() - ! Collect results and statistics call calculate_keff() - ! print cycle information - ! Turn tallies on once inactive cycles are complete - if (current_cycle == n_inactive) then + if (current_batch == n_inactive) then tallies_on = .true. call timer_stop(time_inactive) call timer_start(time_active) end if - ! Stop timer for inter-cycle synchronization - call timer_stop(time_intercycle) - - end do CYCLE_LOOP + end do BATCH_LOOP call timer_stop(time_active) diff --git a/src/output.F90 b/src/output.F90 index 1a9141c0e0..920d9643f7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -830,8 +830,9 @@ contains call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) if (problem_type == PROB_CRITICALITY) then write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' - write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles - write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive + write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches + write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive + write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch elseif (problem_type == PROB_SOURCE) then write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' end if @@ -950,9 +951,9 @@ contains write(ou,100) "Total time in simulation", time_inactive % elapsed + & time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed - write(ou,100) " Time in inactive cycles", time_inactive % elapsed - write(ou,100) " Time in active cycles", time_active % elapsed - write(ou,100) " Time between cycles", time_intercycle % elapsed + write(ou,100) " Time in inactive batches", time_inactive % elapsed + write(ou,100) " Time in active batches", time_active % elapsed + write(ou,100) " Time between generations", time_intercycle % elapsed write(ou,100) " Accumulating tallies", time_ic_tallies % elapsed write(ou,100) " Sampling source sites", time_ic_sample % elapsed write(ou,100) " SEND/RECV source sites", time_ic_sendrecv % elapsed @@ -960,7 +961,7 @@ contains write(ou,100) "Total time elapsed", time_total % elapsed ! display calculate rate and final keff - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) string = to_str(speed) diff --git a/src/source.F90 b/src/source.F90 index fd9cbf1c68..943f18c52c 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -145,12 +145,13 @@ contains p % id = bank_first + index_source - 1 ! set random number seed - particle_seed = (current_cycle - 1)*n_particles + p % id + particle_seed = ((current_batch - 1)*gen_per_batch + & + current_gen - 1)*n_particles + p % id call set_particle_seed(particle_seed) ! set particle trace trace = .false. - if (current_cycle == trace_cycle .and. & + if (current_batch == trace_batch .and. current_gen == trace_gen .and. & p % id == trace_particle) trace = .true. end subroutine get_source_particle diff --git a/src/tally.F90 b/src/tally.F90 index af6d93e738..555e873ad2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1048,28 +1048,31 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_CYCLE_ESTIMATE +! ACCUMULATE_BATCH_ESTIMATE !=============================================================================== - elemental subroutine accumulate_cycle_estimate(score) + elemental subroutine accumulate_batch_estimate(score) type(TallyScore), intent(inout) :: score + real(8) :: val + ! Add the sum and square of the sum of contributions from each cycle ! within a cycle to the variables sum and sum_sq. This will later allow us ! to calculate a variance on the tallies - score % sum = score % sum + score % value/n_particles - score % sum_sq = score % sum_sq + (score % value/n_particles)**2 + val = score % value/(n_particles*gen_per_batch) + score % sum = score % sum + val + score % sum_sq = score % sum_sq + val*val - ! Reset the single cycle estimate + ! Reset the single batch estimate score % value = ZERO - end subroutine accumulate_cycle_estimate + end subroutine accumulate_batch_estimate !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history -! within the cycle to a new random variable +! within the batch to a new random variable !=============================================================================== subroutine synchronize_tallies() @@ -1086,7 +1089,7 @@ contains t => tallies(i) ! Loop over all filter and scoring bins - call accumulate_cycle_estimate(t % scores) + call accumulate_batch_estimate(t % scores) end do end subroutine synchronize_tallies diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 23247fe991..69388420ca 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -4,9 +4,10 @@ - + + From 47364fae97a9951a03e62fb3f59988bb35d99d70 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Mar 2012 16:53:56 -0400 Subject: [PATCH 12/42] Added option to use non-unionized energy grid. --- src/constants.F90 | 6 ++++++ src/cross_section.F90 | 27 +++++++++++++++++++----- src/global.F90 | 3 ++- src/initialize.F90 | 8 ++++--- src/input_xml.F90 | 15 +++++++++++++ src/xml-fortran/templates/settings_t.xml | 11 +++++----- 6 files changed, 56 insertions(+), 14 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 89490eb183..67e7877467 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -312,6 +312,12 @@ module constants integer, parameter :: ERROR_INT = -huge(0) real(8), parameter :: ERROR_REAL = -huge(0.0_8) * 0.917826354_8 + ! Energy grid methods + integer, parameter :: & + GRID_NUCLIDE = 1, & ! non-unionized energy grid (MCNP) + GRID_UNION = 2, & ! union grid with pointers + GRID_LETHARGY = 3 ! lethargy mapping (MC21) + ! Source types integer, parameter :: & SRC_BOX = 1, & ! Source in a rectangular prism diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 966ff1cd94..b43283e037 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -37,7 +37,7 @@ contains mat => materials(p % material) ! Find energy index on unionized grid - call find_energy_index() + if (grid_method == GRID_UNION) call find_energy_index() ! Check if there's an S(a,b) table for this material if (mat % has_sab_table) then @@ -107,11 +107,28 @@ contains ! Set pointer to nuclide nuc => nuclides(index_nuclide) - ! TODO: If not using unionized energy grid, we need to find the index on the - ! nuclide energy grid using lethargy mapping or whatever other technique + ! Determine index on nuclide energy grid + select case (grid_method) + case (GRID_UNION) + ! If we're using the unionized grid with pointers, finding the index on + ! the nuclide energy grid is as simple as looking up the pointer - ! get index on nuclide energy grid - IE = nuc % grid_index(p % IE) + IE = nuc % grid_index(p % IE) + + case (GRID_NUCLIDE) + ! If we're not using the unionized grid, we have to do a binary search on + ! the nuclide energy grid in order to determine which points to + ! interpolate between + + if (p % E < nuc % energy(1)) then + IE = 1 + elseif (p % E > nuc % energy(nuc % n_grid)) then + IE = nuc % n_grid - 1 + else + IE = binary_search(nuc % energy, nuc % n_grid, p % E) + end if + + end select ! check for rare case where two energy points are the same if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 diff --git a/src/global.F90 b/src/global.F90 index d7edf18215..c64ebedc9d 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -82,7 +82,8 @@ module global type(DictionaryCI), pointer :: xs_listing_dict => null() ! Unionized energy grid - integer :: n_grid ! number of points on unionized grid + integer :: grid_method ! how to treat the energy grid + integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid ! Unreoslved resonance probablity tables diff --git a/src/initialize.F90 b/src/initialize.F90 index 697273533e..3d3719b48d 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -106,9 +106,11 @@ contains call timer_stop(time_read_xs) ! Construct unionized energy grid from cross-sections - call timer_start(time_unionize) - call unionized_grid() - call timer_stop(time_unionize) + if (grid_method == GRID_UNION) then + call timer_start(time_unionize) + call unionized_grid() + call timer_stop(time_unionize) + end if ! Create tally map call create_tally_map() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 29c72309a9..c1257855e1 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -67,6 +67,7 @@ contains ! Initialize XML scalar variables cross_sections_ = "" verbosity_ = 0 + energy_grid_ = "union" ! Parse settings.xml file call read_xml_file_settings_t(filename) @@ -111,6 +112,20 @@ contains call fatal_error() end if + ! Energy grid methods + select case (energy_grid_) + case ('nuclide') + grid_method = GRID_NUCLIDE + case ('union') + grid_method = GRID_UNION + case ('lethargy') + message = "Lethargy mapped energy grid not yet supported." + call fatal_error() + case default + message = "Unknown energy grid method: " // energy_grid_ + call fatal_error() + end select + ! Verbosity if (verbosity_ > 0) verbosity = verbosity_ diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 69388420ca..4b450f5724 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -27,13 +27,14 @@ - + + + + + - - - + - From ad53071fc417085cfa3cc3639a03c3c7b8ccffc7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 10:43:27 -0400 Subject: [PATCH 13/42] Updated documentation for batches and energy_grid. --- docs/source/usersguide/input.rst | 31 +++++++++++++++++++++++++------ 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index e213d65c48..55072f24fa 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -307,16 +307,23 @@ settings.xml file. The ``criticality`` element indicates that a criticality calculation should be performed. It has the following attributes/sub-elements: - :cycles: - The number of total fission source iterations. + :batches: + The total number of batches, where each batch corresponds to multiple + fission source iterations. Batching is done to eliminate correlation between + realizations of random variables. *Default*: None + :generations_per_batch: + The number of total fission source iterations per batch. + + *Default*: 1 + :inactive: - The number of inactive fission source iterations. In general, the starting - cycles in a criticality calculation can not be used to contribute to tallies - since the fission source distribution and eigenvalue are generally not - converged immediately + The number of inactive batches. In general, the starting cycles in a + criticality calculation can not be used to contribute to tallies since the + fission source distribution and eigenvalue are generally not converged + immediately. *Default*: None @@ -353,6 +360,18 @@ default. This element has the following attributes/sub-elements: *Default*: 1.0 +``energy_grid`` Element +----------------------- + +The ``energy_grid`` element determines the treatment of the energy grid during a +simulation. Setting this element to "nuclide" will cause OpenMC to use a +nuclide's energy grid when determining what points to interpolate between for +determining cross sections (i.e. non-unionized energy grid). To use a unionized +energy grid, set this element to "union". Note that the unionized energy grid +treatment is slightly different than that employed in Serpent. + + *Default*: union + ``entropy`` Element ------------------- From 09f6ea64034c497649256792bfc8bf5b46715d87 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 10:44:42 -0400 Subject: [PATCH 14/42] Updated examples to use batches instead of cycles. --- examples/basic/settings.xml | 2 +- examples/lattice/nested/settings.xml | 2 +- examples/lattice/simple/settings.xml | 2 +- examples/reflective/settings.xml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/basic/settings.xml b/examples/basic/settings.xml index e9c2e78429..73dd3cfbe0 100644 --- a/examples/basic/settings.xml +++ b/examples/basic/settings.xml @@ -3,7 +3,7 @@ - 15 + 15 5 10000 diff --git a/examples/lattice/nested/settings.xml b/examples/lattice/nested/settings.xml index aa0b98ee4a..2c440b3049 100644 --- a/examples/lattice/nested/settings.xml +++ b/examples/lattice/nested/settings.xml @@ -3,7 +3,7 @@ - 20 + 20 10 10000 diff --git a/examples/lattice/simple/settings.xml b/examples/lattice/simple/settings.xml index aa0b98ee4a..2c440b3049 100644 --- a/examples/lattice/simple/settings.xml +++ b/examples/lattice/simple/settings.xml @@ -3,7 +3,7 @@ - 20 + 20 10 10000 diff --git a/examples/reflective/settings.xml b/examples/reflective/settings.xml index 02ad736305..978cb25ed9 100644 --- a/examples/reflective/settings.xml +++ b/examples/reflective/settings.xml @@ -3,7 +3,7 @@ - 500 + 500 10 1000 From c2e6fda374643949f2a3a5af7fc3a9741ad531d4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 12:46:53 -0400 Subject: [PATCH 15/42] Incremented version number to 0.4.1. --- src/constants.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index 67e7877467..c54d7afe44 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -7,7 +7,7 @@ module constants integer, parameter :: VERSION_MAJOR = 0 integer, parameter :: VERSION_MINOR = 4 - integer, parameter :: VERSION_RELEASE = 0 + integer, parameter :: VERSION_RELEASE = 1 ! ============================================================================ ! ADJUSTABLE PARAMETERS From 516f7008f697ba0ce33e20b8a328492227add604 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 13:36:47 -0400 Subject: [PATCH 16/42] Added release notes for versions 0.4.0 and 0.4.1. --- docs/source/index.rst | 1 + docs/source/releasenotes/index.rst | 14 ++++++ docs/source/releasenotes/notes_0.4.0.rst | 36 ++++++++++++++++ docs/source/releasenotes/notes_0.4.1.rst | 55 ++++++++++++++++++++++++ 4 files changed, 106 insertions(+) create mode 100644 docs/source/releasenotes/index.rst create mode 100644 docs/source/releasenotes/notes_0.4.0.rst create mode 100644 docs/source/releasenotes/notes_0.4.1.rst diff --git a/docs/source/index.rst b/docs/source/index.rst index 828f5335d1..18db6e8628 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -25,6 +25,7 @@ Contents :maxdepth: 1 install + releasenotes/index methods/index usersguide/index developers diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst new file mode 100644 index 0000000000..419c1dc4d9 --- /dev/null +++ b/docs/source/releasenotes/index.rst @@ -0,0 +1,14 @@ +.. _releasenotes: + +============= +Release Notes +============= + +The release notes for OpenMC give a list of system requirements, new features, +bugs fixed, and known issues for each successive release. + +.. toctree:: + :maxdepth: 1 + + notes_0.4.1 + notes_0.4.0 diff --git a/docs/source/releasenotes/notes_0.4.0.rst b/docs/source/releasenotes/notes_0.4.0.rst new file mode 100644 index 0000000000..3e0a444e35 --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.0.rst @@ -0,0 +1,36 @@ +.. _notes_0.4.0: + +============================== +Release Notes for OpenMC 0.4.0 +============================== + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions as well as +Mac OS X. However, it has not been tested yet on any releases of Microsoft +Windows. Memory requirements will vary depending on the size of the problem at +hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- The probability table method for treatment of energy self-shielding in the + unresolved resonance range has been implemented and is now turned on by + default. +- Calculation of Shannon entropy for assessing convergence of the fission source + distribution. +- Ability to compile with the PGI Fortran compiler. +- Ability to run on IBM BlueGene/P machines. +- Completely rewrote how nested universes are handled. Geometry is now much more + robust. + +--------- +Bug Fixes +--------- + +- Many geometry errors have been fixed. The Monte Carlo performance benchmark + can now be successfully run in OpenMC. diff --git a/docs/source/releasenotes/notes_0.4.1.rst b/docs/source/releasenotes/notes_0.4.1.rst new file mode 100644 index 0000000000..2492f23a98 --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.1.rst @@ -0,0 +1,55 @@ +.. _notes_0.4.1: + +============================== +Release Notes for OpenMC 0.4.1 +============================== + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions as well as +Mac OS X. However, it has not been tested yet on any releases of Microsoft +Windows. Memory requirements will vary depending on the size of the problem at +hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- A batching method has been implemented so that statistics can be calculated + based on multiple generations instead of a single generation. This can help to + overcome problems with underpredicted variance in problems where there is + correlation between successive fission source iterations. +- Users now have the option to select a non-unionized energy grid for problems + with many nuclides where the use of a unionized grid is not feasible. +- Improved plotting capability (Nick Horelik). The plotting input is now in + ``plots.xml`` instead of ``plot.xml``. +- Added multiple estimators for k-effective and added a global tally for + leakage. +- Moved cross section-related output into cross_sections.out. +- Improved timing capabilities. +- Can now use more than 2**31 - 1 particles per generation. +- Improved fission bank synchronization method. This also necessitated changing + the source bank to be of type Bank rather than of type Particle. +- Added HDF5 output (not complete yet). +- Major changes to tally implementation. + +--------- +Bug Fixes +--------- + +- `b206a8`_: Fixed subtle error in the sampling of energy distributions. +- `800742`_: Fixed error in sampling of angle and rotating angles. +- `a07c08`_: Fixed bug in linear-linear interpolation during sampling energy. +- `a75283`_: Fixed energy and energyout tally filters to support many bins. +- `95cfac`_: Fixed error in cell neighbor searches. +- `83a803`_: Fixed bug related to probability tables. + +.. _b206a8: https://github.com/mit-crpg/openmc/commit/b206a8 +.. _800742: https://github.com/mit-crpg/openmc/commit/800742 +.. _a07c08: https://github.com/mit-crpg/openmc/commit/a07c08 +.. _a75283: https://github.com/mit-crpg/openmc/commit/a75283 +.. _95cfac: https://github.com/mit-crpg/openmc/commit/95cfac +.. _83a803: https://github.com/mit-crpg/openmc/commit/83a803 From 838bb351fcf7ac6b75287475eedb898f508a8843 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 14:07:27 -0400 Subject: [PATCH 17/42] Added option to change random number seed based on cmfd_parallel branch. --- docs/source/usersguide/input.rst | 8 ++++++++ src/global.F90 | 3 +++ src/initialize.F90 | 6 +++--- src/input_xml.F90 | 4 ++++ src/random_lcg.F90 | 4 +++- src/xml-fortran/templates/settings_t.xml | 1 + 6 files changed, 22 insertions(+), 4 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 55072f24fa..3e0aec34e7 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -404,6 +404,14 @@ sub-elements and can be set to either "off" or "on". *Default*: on +``seed`` Element +---------------- + +The ``seed`` element is used to set the seed used for the linear congruential +pseudo-random number generator. + + *Default*: 1 + ``source`` Element ------------------ diff --git a/src/global.F90 b/src/global.F90 index c64ebedc9d..5fadb2e1af 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -211,6 +211,9 @@ module global ! Message used in message/warning/fatal_error character(MAX_LINE_LEN) :: message + ! Random number seed + integer(8) :: seed = 1_8 + ! Problem type integer :: problem_type = PROB_CRITICALITY diff --git a/src/initialize.F90 b/src/initialize.F90 index 3d3719b48d..3275c1c282 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -72,15 +72,15 @@ contains call header("INITIALIZATION", level=1) end if - ! Initialize random number generator - call initialize_prng() - ! set up dictionaries call create_dictionaries() ! Read XML input files call read_input_xml() + ! Initialize random number generator + call initialize_prng() + ! Set up universe structures call prepare_universes() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c1257855e1..2837fd1ebf 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -68,6 +68,7 @@ contains cross_sections_ = "" verbosity_ = 0 energy_grid_ = "union" + seed_ = 0_8 ! Parse settings.xml file call read_xml_file_settings_t(filename) @@ -112,6 +113,9 @@ contains call fatal_error() end if + ! Copy random number seed if specified + if (seed_ > 0) seed = seed_ + ! Energy grid methods select case (energy_grid_) case ('nuclide') diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 674be52122..694885c013 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -49,7 +49,9 @@ contains subroutine initialize_prng() - prn_seed0 = 1_8 + use global, only: seed + + prn_seed0 = seed prn_seed = prn_seed prn_mult = 2806196910506780709_8 prn_add = 1_8 diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 4b450f5724..dbe8bd6226 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -32,6 +32,7 @@ + From c9e62980a09ae86bf9cc9f88fd0190083038deec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 14:55:05 -0400 Subject: [PATCH 18/42] Added count_fission_sites in mesh module and restructured Shannon entropy calculation. --- src/DEPENDENCIES | 3 ++ src/global.F90 | 4 +-- src/intercycle.F90 | 72 ++++++++++++----------------------------- src/mesh.F90 | 80 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 105 insertions(+), 54 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index f2cb4f5e12..587fb06cca 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -148,7 +148,10 @@ main.o: string.o main.o: tally.o main.o: timing.o +mesh.o: constants.o +mesh.o: global.o mesh.o: mesh_header.o +mesh.o: particle_header.o output.o: constants.o output.o: datatypes.o diff --git a/src/global.F90 b/src/global.F90 index 5fadb2e1af..752e0db3a7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -152,8 +152,8 @@ module global ! Shannon entropy logical :: entropy_on = .false. - real(8) :: entropy ! value of shannon entropy - real(8), allocatable :: entropy_p(:) ! fraction of source sites in each cell + real(8) :: entropy ! value of shannon entropy + real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh ! ============================================================================ diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 6d3ae131e8..97a7b9f86d 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -4,7 +4,7 @@ module intercycle use error, only: fatal_error, warning use global - use mesh, only: get_mesh_bin + use mesh, only: count_fission_sites use mesh_header, only: StructuredMesh use output, only: write_message use random_lcg, only: prn, set_particle_seed, prn_skip @@ -283,12 +283,10 @@ contains subroutine shannon_entropy() - integer :: i ! index for bank sites + integer :: i, j, k ! index for bank sites integer :: n ! # of boxes in each dimension - integer :: bin ! index in entropy_p - integer(8) :: total_bank ! total # of fission bank sites - integer, save :: n_box ! total # of boxes on mesh - logical :: outside_box ! were there sites outside entropy box? + real(8) :: total ! total weight of fission bank sites + logical :: sites_outside ! were there sites outside entropy box? type(StructuredMesh), pointer :: m => null() ! Get pointer to entropy mesh @@ -312,69 +310,39 @@ contains m % dimension = n end if - ! Determine total number of mesh boxes - n_box = product(m % dimension) - ! allocate and determine width allocate(m % width(3)) m % width = (m % upper_right - m % lower_left) / m % dimension ! allocate p - allocate(entropy_p(n_box)) + allocate(entropy_p(m % dimension(1), m % dimension(2), & + m % dimension(3))) end if - ! initialize p - entropy_p = ZERO - outside_box = .false. - - ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do i = 1, int(n_bank,4) - ! determine scoring bin for entropy mesh - call get_mesh_bin(m, fission_bank(i) % xyz, bin) - - ! if outside mesh, skip particle - if (bin == NO_BIN_FOUND) then - outside_box = .true. - cycle - end if - - ! add to appropriate mesh box - entropy_p(bin) = entropy_p(bin) + 1 - end do FISSION_SITES + ! count number of fission sites over mesh + call count_fission_sites(m, entropy_p, total, sites_outside) ! display warning message if there were sites outside entropy box - if (outside_box) then + if (sites_outside) then message = "Fission source site(s) outside of entropy box." call warning() end if -#ifdef MPI - ! collect values from all processors - if (master) then - call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - else - call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) - end if - - ! determine total number of bank sites - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - ! sum values to obtain shannon entropy if (master) then - ! Normalize to number of bank sites - entropy_p = entropy_p / total_bank + ! Normalize to total weight of bank sites + entropy_p = entropy_p / total entropy = 0 - do i = 1, n_box - if (entropy_p(i) > 0) then - entropy = entropy - entropy_p(i) * log(entropy_p(i))/log(2.0) - end if + do i = 1, m % dimension(1) + do j = 1, m % dimension(2) + do k = 1, m % dimension(3) + if (entropy_p(i,j,k) > ZERO) then + entropy = entropy - entropy_p(i,j,k) * & + log(entropy_p(i,j,k))/log(2.0) + end if + end do + end do end do end if diff --git a/src/mesh.F90 b/src/mesh.F90 index 22c2f28c86..756525327c 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,9 +1,14 @@ module mesh use constants + use global use mesh_header use particle_header, only: Particle +#ifdef MPI + use mpi +#endif + implicit none contains @@ -137,4 +142,79 @@ contains end subroutine bin_to_mesh_indices +!=============================================================================== +! COUNT_FISSION_SITES determines the number of fission bank sites in each cell +! of a given mesh. This can be used for a variety of purposes (Shannon entropy, +! CMFD, uniform fission source weighting) +!=============================================================================== + + subroutine count_fission_sites(m, count, total, sites_outside) + + type(StructuredMesh), pointer :: m ! mesh to count sites + real(8), intent(out) :: count(:,:,:) ! weight of sites in each cell + real(8), intent(out) :: total ! total weight of sites + logical, optional :: sites_outside ! were there sites outside mesh? + + integer :: i ! loop index for local fission sites + integer :: ijk(3) ! indices on mesh + real(8) :: weight ! accumulated weight of sites + logical :: in_mesh ! was single site outside mesh? + logical :: outside ! was any site outside mesh? +#ifdef MPI + integer :: n ! total size of count variable +#endif + + ! initialize variables + count = ZERO + weight = ZERO + outside = .false. + + ! loop over fission sites and count how many are in each mesh box + FISSION_SITES: do i = 1, int(n_bank,4) + ! determine scoring bin for entropy mesh + call get_mesh_indices(m, fission_bank(i) % xyz, ijk, in_mesh) + + ! if outside mesh, skip particle + if (.not. in_mesh) then + outside = .true. + cycle + end if + + ! add weight + weight = weight + ONE + + ! add to appropriate mesh box + ! TODO: if tracking weight through bank, add weight instead + count(ijk(1),ijk(2),ijk(3)) = count(ijk(1),ijk(2),ijk(3)) + 1 + end do FISSION_SITES + +#ifdef MPI + ! determine total number of mesh cells + n = size(count,1) * size(count,2) * size(count,3) + + ! collect values from all processors + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(count, count, n, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! Check if there were sites outside the mesh for any processor + if (present(sites_outside)) then + call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & + MPI_COMM_WORLD, mpi_err) + end if + + ! determine total weight of bank sites + call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total = weight + sites_outside = outside +#endif + + end subroutine count_fission_sites + end module mesh From 2a5b9cdcd518baa1c6cd9b64a3b65d28fd5b4266 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 20 Mar 2012 14:49:14 -0400 Subject: [PATCH 19/42] Fixed regression in plotting. --- src/initialize.F90 | 6 ++- src/input_xml.F90 | 97 +++++++++++++++++++++++----------------------- 2 files changed, 54 insertions(+), 49 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 3275c1c282..6bc2e899c8 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,7 +12,7 @@ module initialize use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE use global use input_xml, only: read_input_xml, read_cross_sections_xml, & - cells_in_univ_dict + cells_in_univ_dict, read_plots_xml use output, only: title, header, print_summary, print_geometry, & print_plot, create_summary_file, & create_xs_summary_file @@ -81,6 +81,10 @@ contains ! Initialize random number generator call initialize_prng() + ! 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 (plotting) call read_plots_xml() + ! Set up universe structures call prepare_universes() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2837fd1ebf..cbfef818a3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -32,7 +32,6 @@ contains call read_geometry_xml() call read_materials_xml() call read_tallies_xml() - if (plotting) call read_plots_xml() end subroutine read_input_xml @@ -1163,14 +1162,16 @@ contains end if ! Copy plot background color - if (size(plot_(i) % background) == 3) then - pl % not_found % rgb = plot_(i) % background - else if (size(plot_(i) % background) == 0) then - pl % not_found % rgb = (/ 255, 255, 255 /) + if (associated(plot_(i) % background)) then + if (size(plot_(i) % background) == 3) then + pl % not_found % rgb = plot_(i) % background + else + message = "Bad background RGB " & + // "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if else - message = "Bad background RGB " & - // "in plot " // trim(to_str(pl % id)) - call fatal_error() + pl % not_found % rgb = (/ 255, 255, 255 /) end if ! Copy plot type @@ -1232,7 +1233,7 @@ contains pl % colors(j) % rgb(3) = prn()*255 end do - case ("mat") + case ("mat", "material") pl % color_by = PLOT_COLOR_MATS allocate(pl % colors(n_materials)) @@ -1249,55 +1250,55 @@ contains end select ! Copy user specified colors - n_cols = size(plot_(i) % col_spec_) - do j = 1, n_cols - if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then - message = "Bad RGB " & + if (associated(plot_(i) % col_spec_)) then + n_cols = size(plot_(i) % col_spec_) + do j = 1, n_cols + if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then + message = "Bad RGB " & // "in plot " // trim(to_str(pl % id)) - call fatal_error() - end if + call fatal_error() + end if - col_id = plot_(i) % col_spec_(j) % id + col_id = plot_(i) % col_spec_(j) % id - if (pl % color_by == PLOT_COLOR_CELLS) then + if (pl % color_by == PLOT_COLOR_CELLS) then - if (dict_has_key(cell_dict, col_id)) then - pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb - else - message = "Could not find cell " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) - call fatal_error() - end if + if (dict_has_key(cell_dict, col_id)) then + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb + else + message = "Could not find cell " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if - else if (pl % color_by == PLOT_COLOR_MATS) then + else if (pl % color_by == PLOT_COLOR_MATS) then - if (dict_has_key(material_dict, col_id)) then - pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb - else - message = "Could not find material " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) - call fatal_error() - end if + if (dict_has_key(material_dict, col_id)) then + pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb + else + message = "Could not find material " // trim(to_str(col_id)) // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + end if - end if - - end do + end if + end do + end if ! Alter colors based on mask information - if (size(plot_(i) % mask_) > 1) then - message = "Mutliple masks" // & - " specified in plot " // trim(to_str(pl % id)) - call fatal_error() - else if (.not. size(plot_(i) % mask_) == 0) then - do j=1,size(pl % colors) - if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then - pl % colors(j) % rgb = plot_(i) % mask_(1) % background - end if - end do + if (associated(plot_(i) % mask_)) then + if (size(plot_(i) % mask_) > 1) then + message = "Mutliple masks" // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() + else if (.not. size(plot_(i) % mask_) == 0) then + do j=1,size(pl % colors) + if (.not. any(j .eq. plot_(i) % mask_(1) % components)) then + pl % colors(j) % rgb = plot_(i) % mask_(1) % background + end if + end do + end if end if - - - end do From 24bed2624782872e654f0baeeffc9c5f925b02d1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Mar 2012 17:45:33 -0400 Subject: [PATCH 20/42] Added point type for starting source. --- src/constants.F90 | 3 +-- src/input_xml.F90 | 3 +++ src/source.F90 | 2 ++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index c54d7afe44..be0f712fc7 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -321,8 +321,7 @@ module constants ! Source types integer, parameter :: & SRC_BOX = 1, & ! Source in a rectangular prism - SRC_CELL = 2, & ! Source in a cell - SRC_SURFACE = 3 ! Source on a surface + SRC_POINT = 2 ! Source at a single point integer, parameter :: & PROB_SOURCE = 1, & ! External source problem diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cbfef818a3..3f92b0b250 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -140,6 +140,9 @@ contains case ('box') external_source % type = SRC_BOX coeffs_reqd = 6 + case ('point') + external_source % type = SRC_POINT + coeffs_reqd = 3 case default message = "Invalid source type: " // trim(source_ % type) call fatal_error() diff --git a/src/source.F90 b/src/source.F90 index 943f18c52c..e9088e3ddb 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -93,6 +93,8 @@ contains p_max = external_source % values(4:6) r = (/ (prn(), k = 1,3) /) source_bank(j) % xyz = p_min + r*(p_max - p_min) + case (SRC_POINT) + source_bank(j) % xyz = external_source % values end select ! sample angle From c8b37112cd0724a40dcd1a83920df010d7985d9f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Mar 2012 11:41:45 -0400 Subject: [PATCH 21/42] Added units='sum' option for specifying density. --- src/input_xml.F90 | 65 +++++++++++++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3f92b0b250..b5e29a5331 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -518,6 +518,7 @@ contains integer :: n ! number of nuclides real(8) :: val ! value entered for density logical :: file_exists ! does materials.xml exist? + logical :: sum_density ! density is taken to be sum of nuclide densities character(3) :: default_xs ! default xs identifier (e.g. 70c) character(12) :: name ! name of nuclide character(MAX_WORD_LEN) :: units ! units on density @@ -557,32 +558,43 @@ contains ! Copy material id m % id = material_(i) % id - ! Copy density -- the default value for the units is given in the - ! material_t.xml file and doesn't need to be specified here, hence case - ! default results in an error. - val = material_(i) % density % value - if (val <= ZERO) then - message = "Need to specify a positive density on material " // & - trim(to_str(m % id)) // "." - call fatal_error() - end if - + ! Copy value and units + val = material_(i) % density % value units = material_(i) % density % units - call lower_case(units) - select case(trim(units)) - case ('g/cc', 'g/cm3') - m % density = -val - case ('kg/m3') - m % density = -0.001 * val - case ('atom/b-cm') - m % density = val - case ('atom/cm3', 'atom/cc') - m % density = 1.0e-24 * val - case default - message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(to_str(m % id)) - call fatal_error() - end select + + if (units == 'sum') then + ! If the user gave the units as 'sum', then the total density of the + ! material is taken to be the sum of the atom fractions listed on the + ! nuclides + + sum_density = .true. + + else + ! Check for erroneous density + sum_density = .false. + if (val <= ZERO) then + message = "Need to specify a positive density on material " // & + trim(to_str(m % id)) // "." + call fatal_error() + end if + + ! Adjust material density based on specified units + call lower_case(units) + select case(trim(units)) + case ('g/cc', 'g/cm3') + m % density = -val + case ('kg/m3') + m % density = -0.001 * val + case ('atom/b-cm') + m % density = val + case ('atom/cm3', 'atom/cc') + m % density = 1.0e-24 * val + case default + message = "Unkwown units '" // trim(material_(i) % density % units) & + // "' specified on material " // trim(to_str(m % id)) + call fatal_error() + end select + end if ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then @@ -657,6 +669,9 @@ contains end if end do + ! Determine density if it is a sum value + if (sum_density) m % density = sum(m % atom_percent) + ! Add material to dictionary call dict_add_key(material_dict, m % id, i) From 20aac87d607d85abe16417d734750ee2ad134fbb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Mar 2012 12:19:41 -0400 Subject: [PATCH 22/42] Added preliminary release notes for 0.4.2. --- docs/source/releasenotes/index.rst | 1 + docs/source/releasenotes/notes_0.4.2.rst | 37 ++++++++++++++++++++++++ 2 files changed, 38 insertions(+) create mode 100644 docs/source/releasenotes/notes_0.4.2.rst diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst index 419c1dc4d9..8f6c633421 100644 --- a/docs/source/releasenotes/index.rst +++ b/docs/source/releasenotes/index.rst @@ -10,5 +10,6 @@ bugs fixed, and known issues for each successive release. .. toctree:: :maxdepth: 1 + notes_0.4.2 notes_0.4.1 notes_0.4.0 diff --git a/docs/source/releasenotes/notes_0.4.2.rst b/docs/source/releasenotes/notes_0.4.2.rst new file mode 100644 index 0000000000..8d4c1f2d2a --- /dev/null +++ b/docs/source/releasenotes/notes_0.4.2.rst @@ -0,0 +1,37 @@ +.. _notes_0.4.2: + +============================== +Release Notes for OpenMC 0.4.2 +============================== + +.. note:: + These release notes are for an upcoming release of OpenMC and are still + subject to change. + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions as well as +Mac OS X. However, it has not been tested yet on any versions of Microsoft +Windows. Memory requirements will vary depending on the size of the problem at +hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- Ability to specify 'point' source type. +- Ability to change random number seed. +- User's can now specify units='sum' on a tag. This tells the code + that the total material density is the sum of the atom fractions listed for + each nuclide on the material. + +--------- +Bug Fixes +--------- + +- `2a5b9c`_: Fixed regression in plotting. + +.. _2a5b9c: https://github.com/mit-crpg/openmc/commit/2a5b9c From d050c72d1d7e12df6fdc4d2ede24f481950f2265 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Mar 2012 14:16:29 -0400 Subject: [PATCH 23/42] Added Bessel's correction to make estimate of variance unbiased. --- src/intercycle.F90 | 7 +++++-- src/tally.F90 | 8 ++++++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 97a7b9f86d..3d5f0cd704 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -394,12 +394,15 @@ contains ! Accumulate single batch realizations of k call accumulate_batch_estimate(global_tallies) - ! Determine mean and standard deviation of mean + ! Determine sample mean of keff keff = global_tallies(K_ANALOG) % sum/n - keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n) ! Display output for this batch if (current_batch > n_inactive + 1) then + ! Determine standard deviation of the sample mean of keff + keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - & + keff*keff)/(n - 1)) + if (entropy_on) then write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std diff --git a/src/tally.F90 b/src/tally.F90 index 555e873ad2..4702ed9baf 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1532,9 +1532,13 @@ contains type(TallyScore), intent(inout) :: score - ! Calculate mean and standard deviation of the mean + ! Calculate sample mean and standard deviation of the mean -- note that we + ! have used Bessel's correction so that the estimator of the variance of the + ! sample mean is unbiased. + score % sum = score % sum/n_active - score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2)/n_active) + score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & + (n_active - 1)) end subroutine calculate_statistics From 28aa9111de1776adf393121d8c53af98a9182e2b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 10:16:13 -0400 Subject: [PATCH 24/42] Added trace message for cells. --- src/geometry.F90 | 6 ++++++ src/physics.F90 | 16 +++++----------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 3457719556..fe75d2c4e6 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -136,6 +136,12 @@ contains ! Set cell on this level p % coord % cell = index_cell + ! Show cell information on trace + if (trace) then + message = " Entering cell " // trim(to_str(c % id)) + call write_message() + end if + if (c % type == CELL_NORMAL) then ! ================================================================= ! AT LOWEST UNIVERSE, TERMINATE SEARCH diff --git a/src/physics.F90 b/src/physics.F90 index 8ed65c8c7f..a553a115d8 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -39,6 +39,11 @@ contains logical :: found_cell ! found cell which particle is in? type(LocalCoord), pointer :: coord => null() + if (verbosity >= 9 .or. trace) then + message = "Simulating Particle " // trim(to_str(p % id)) + call write_message() + end if + if (p % coord % cell == NONE) then call find_cell(found_cell) ! Particle couldn't be located @@ -51,17 +56,6 @@ contains p % cell_born = p % coord % cell end if - if (verbosity >= 9 .or. trace) then - message = "Simulating Particle " // trim(to_str(p % id)) - call write_message() - end if - - if (verbosity >= 10 .or. trace) then - message = " Born in cell " // trim(to_str(& - cells(p % coord % cell) % id)) - call write_message() - end if - ! Initialize number of events to zero n_event = 0 From a272bafb03d0c4ea1b5dee5bb41649e99691da42 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 10:38:11 -0400 Subject: [PATCH 25/42] Added trace information on collision. --- src/geometry.F90 | 2 +- src/particle_header.F90 | 3 ++- src/physics.F90 | 18 +++++++++++++++--- 3 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index fe75d2c4e6..0abda64c13 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -137,7 +137,7 @@ contains p % coord % cell = index_cell ! Show cell information on trace - if (trace) then + if (verbosity >= 10 .or. trace) then message = " Entering cell " // trim(to_str(c % id)) call write_message() end if diff --git a/src/particle_header.F90 b/src/particle_header.F90 index a9cfa12c24..4330dfb882 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -52,7 +52,8 @@ module particle_header real(8) :: last_E ! last energy ! What event last took place - integer :: event + integer :: event ! scatter, absorption, fission + integer :: event_nuclide ! index in nuclides array ! Post-collision physical data integer :: n_bank ! number of fission sites banked diff --git a/src/physics.F90 b/src/physics.F90 index a553a115d8..27ec0adc20 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -182,9 +182,18 @@ contains ! Display information about collision if (verbosity >= 10 .or. trace) then -!!$ message = " " // trim(reaction_name(MT)) // ". Energy = " // & -!!$ trim(to_str(p % E * 1e6_8)) // " eV." -!!$ call write_message() + select case (p % event) + case (EVENT_SCATTER) + message = " Scattered off of" + case (EVENT_FISSION) + message = " Fissioned in" + case (EVENT_ABSORB) + message = " Absorbed in" + end select + message = trim(message) // " " // trim(nuclides(& + p % event_nuclide) % name) // ". Energy = " // & + trim(to_str(p % E * 1e6_8)) // " eV." + call write_message() end if ! check for very low energy @@ -262,6 +271,9 @@ contains IE = micro_xs(index_nuclide) % index_grid f = micro_xs(index_nuclide) % interp_factor + ! Save which nuclide particle had collision with + p % event_nuclide = index_nuclide + ! ========================================================================== ! DISAPPEARANCE REACTIONS (ANALOG) OR IMPLICIT CAPTURE (SURVIVAL BIASING) From 057f09444aa49f25f0f301c15fc77604094ba327 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 13:58:48 -0400 Subject: [PATCH 26/42] Added estimator for diffusion coefficient. --- src/constants.F90 | 19 ++++++++++--------- src/input_xml.F90 | 5 +++++ src/tally.F90 | 15 +++++++++++++++ 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index be0f712fc7..9cb5c99d55 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -243,7 +243,7 @@ module constants EVENT_FISSION = 3 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 15 + integer, parameter :: N_SCORE_TYPES = 16 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -252,14 +252,15 @@ module constants SCORE_SCATTER_1 = -5, & ! first scattering moment SCORE_SCATTER_2 = -6, & ! second scattering moment SCORE_SCATTER_3 = -7, & ! third scattering moment - SCORE_N_1N = -8, & ! (n,1n) rate - SCORE_N_2N = -9, & ! (n,2n) rate - SCORE_N_3N = -10, & ! (n,3n) rate - SCORE_N_4N = -11, & ! (n,4n) rate - SCORE_ABSORPTION = -12, & ! absorption rate - SCORE_FISSION = -13, & ! fission rate - SCORE_NU_FISSION = -14, & ! neutron production rate - SCORE_CURRENT = -15 ! partial current + SCORE_DIFFUSION = -8, & ! diffusion coefficient + SCORE_N_1N = -9, & ! (n,1n) rate + SCORE_N_2N = -10, & ! (n,2n) rate + SCORE_N_3N = -11, & ! (n,3n) rate + SCORE_N_4N = -12, & ! (n,4n) rate + SCORE_ABSORPTION = -13, & ! absorption rate + SCORE_FISSION = -14, & ! fission rate + SCORE_NU_FISSION = -15, & ! neutron production rate + SCORE_CURRENT = -16 ! partial current ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b5e29a5331..bb9841a14e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1006,6 +1006,11 @@ contains case ('scatter-3') t % score_bins(j) % scalar = SCORE_SCATTER_3 + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + case ('diffusion') + t % score_bins(j) % scalar = SCORE_DIFFUSION + ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('n1n') diff --git a/src/tally.F90 b/src/tally.F90 index 4702ed9baf..db57a56802 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -345,6 +345,20 @@ contains score = last_wgt * 0.5*(5.0*mu*mu*mu - 3.0*mu) + case (SCORE_DIFFUSION) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle + + ! Temporarily store the scattering cross section + score = material_xs % total - material_xs % absorption + + ! Since this only gets tallied at every scattering event, the flux + ! estimator is 1/Sigma_s. Therefore, the diffusion coefficient + ! times flux is 1/(3*Sigma_s*(Sigma_t - mu*Sigma_s)). + + score = last_wgt / (3.0_8 * score * (material_xs % total - & + mu * score)) + case (SCORE_N_1N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle @@ -1186,6 +1200,7 @@ contains score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment" score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment" score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment" + score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient" score_name(abs(SCORE_N_1N)) = "(n,1n) Rate" score_name(abs(SCORE_N_2N)) = "(n,2n) Rate" score_name(abs(SCORE_N_3N)) = "(n,3n) Rate" From dffed0bfda56bf2d79671effb34776c7dc3d228a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 14:31:12 -0400 Subject: [PATCH 27/42] Added ability to write out source bank at end of simulation (not in parallel yet). --- src/DEPENDENCIES | 1 + src/constants.F90 | 1 + src/finalize.F90 | 6 ++++- src/global.F90 | 3 +++ src/input_xml.F90 | 4 +++ src/source.F90 | 27 ++++++++++++++++++++ src/utils/read_source.py | 32 ++++++++++++++++++++++++ src/xml-fortran/templates/settings_t.xml | 1 + 8 files changed, 74 insertions(+), 1 deletion(-) create mode 100755 src/utils/read_source.py diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 587fb06cca..bcc3b533ff 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -40,6 +40,7 @@ energy_grid.o: output.o finalize.o: global.o finalize.o: hdf5_interface.o finalize.o: output.o +finalize.o: source.o finalize.o: tally.o finalize.o: timing.o diff --git a/src/constants.F90 b/src/constants.F90 index 9cb5c99d55..4edebf1291 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -333,5 +333,6 @@ module constants integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file + integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file end module constants diff --git a/src/finalize.F90 b/src/finalize.F90 index c4dcaf0623..e239c6db9f 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -2,6 +2,7 @@ module finalize use global use output, only: print_runtime + use source, only: write_source_binary use tally, only: write_tallies, tally_statistics use timing, only: timer_start, timer_stop @@ -23,10 +24,13 @@ contains ! Start finalization timer call timer_start(time_finalize) - ! Calculate statistics for tallies and write to tallies.out if (.not. plotting) then + ! Calculate statistics for tallies and write to tallies.out call tally_statistics() if (master) call write_tallies() + + ! Write out binary source + if (write_source) call write_source_binary() end if ! stop timers and show timing statistics diff --git a/src/global.F90 b/src/global.F90 index 752e0db3a7..6cbb99cfa4 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -156,6 +156,9 @@ module global real(8), allocatable :: entropy_p(:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh + ! Write source at end of simulation + logical :: write_source = .false. + ! ============================================================================ ! PARALLEL PROCESSING VARIABLES diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bb9841a14e..f90a54a247 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -68,6 +68,7 @@ contains verbosity_ = 0 energy_grid_ = "union" seed_ = 0_8 + write_source = "" ! Parse settings.xml file call read_xml_file_settings_t(filename) @@ -231,6 +232,9 @@ contains entropy_on = .true. end if + ! Check if the user has specified to write binary source file + if (trim(write_source_) == 'on') write_source = .true. + end subroutine read_settings_xml !=============================================================================== diff --git a/src/source.F90 b/src/source.F90 index e9088e3ddb..0652d22ac3 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -189,4 +189,31 @@ contains end subroutine initialize_particle +!=============================================================================== +! WRITE_SOURCE writes out the final source distribution to a binary file that +! can be used as a starting source in a new simulation +!=============================================================================== + + subroutine write_source_binary() + + integer(8) :: i + + open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', & + ACCESS='stream') + + ! write the number of particles + write(UNIT=UNIT_SOURCE) n_particles + + ! write information from the source bank + do i = 1_8, work + write(UNIT=UNIT_SOURCE) source_bank(i) % xyz + write(UNIT=UNIT_SOURCE) source_bank(i) % uvw + write(UNIT=UNIT_SOURCE) source_bank(i) % E + end do + + ! close binary source file + close(UNIT=UNIT_SOURCE) + + end subroutine write_source_binary + end module source diff --git a/src/utils/read_source.py b/src/utils/read_source.py new file mode 100755 index 0000000000..6da8ff5008 --- /dev/null +++ b/src/utils/read_source.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python + +import struct + +class SourceFile(object): + def __init__(self, filename): + # Open source file for reading + self.f = open(filename, 'r') + + # Read number of source sites + self.n_sites = struct.unpack('l', self.f.read(8))[0] + + # Create list to store source sites + self.sites = [] + + for i in range(self.n_sites): + # Read position, angle, and energy + xyz = struct.unpack('3d', self.f.read(24)) + uvw = struct.unpack('3d', self.f.read(24)) + E = struct.unpack('d', self.f.read(8)) + + # Create source site and append to list + self.sites.append(SourceSite(xyz, uvw, E)) + +class SourceSite(object): + def __init__(self, xyz, uvw, E): + self.xyz = xyz + self.uvw = uvw + self.E = E + + def __repr__(self): + return "".format(self.xyz, self.E) diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index dbe8bd6226..f2b387316d 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -37,5 +37,6 @@ + From 17bc90a52eaa1b6d4b97d07b6aafd4eb271526bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 15:59:01 -0400 Subject: [PATCH 28/42] Ability to start from binary source file (not in parallel). --- src/constants.F90 | 3 ++- src/input_xml.F90 | 26 ++++++++++++++---------- src/source.F90 | 52 ++++++++++++++++++++++++++++++++++++++++------- 3 files changed, 62 insertions(+), 19 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 4edebf1291..90828d3019 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -322,7 +322,8 @@ module constants ! Source types integer, parameter :: & SRC_BOX = 1, & ! Source in a rectangular prism - SRC_POINT = 2 ! Source at a single point + SRC_POINT = 2, & ! Source at a single point + SRC_FILE = 3 ! Source from a file integer, parameter :: & PROB_SOURCE = 1, & ! External source problem diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f90a54a247..d9385fd376 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -137,29 +137,33 @@ contains ! Determine external source type type = source_ % type call lower_case(type) - select case (trim(type)) + select case (type) case ('box') external_source % type = SRC_BOX coeffs_reqd = 6 case ('point') external_source % type = SRC_POINT coeffs_reqd = 3 + case ('file') + external_source % type = SRC_FILE case default message = "Invalid source type: " // trim(source_ % type) call fatal_error() end select ! Coefficients for external surface - n = size(source_ % coeffs) - if (n < coeffs_reqd) then - message = "Not enough coefficients specified for external source." - call fatal_error() - elseif (n > coeffs_reqd) then - message = "Too many coefficients specified for external source." - call fatal_error() - else - allocate(external_source % values(n)) - external_source % values = source_ % coeffs + if (type /= 'file') then + n = size(source_ % coeffs) + if (n < coeffs_reqd) then + message = "Not enough coefficients specified for external source." + call fatal_error() + elseif (n > coeffs_reqd) then + message = "Too many coefficients specified for external source." + call fatal_error() + else + allocate(external_source % values(n)) + external_source % values = source_ % coeffs + end if end if end if diff --git a/src/source.F90 b/src/source.F90 index 0652d22ac3..a9de74bed1 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -43,6 +43,13 @@ contains ! Determine maximum amount of particles to simulate on each processor maxwork = ceiling(real(n_particles)/n_procs,8) + ! ID's of first and last source particles + bank_first = rank*maxwork + 1 + bank_last = min((rank+1)*maxwork, n_particles) + + ! number of particles for this processor + work = bank_last - bank_first + 1 + ! Allocate source bank allocate(source_bank(maxwork), STAT=alloc_err) if (alloc_err /= 0) then @@ -69,16 +76,16 @@ contains call fatal_error() end if + ! Read the source from a binary file instead of sampling from some assumed + ! source distribution + if (external_source % type == SRC_FILE) then + call read_source_binary() + return + end if + ! Initialize first cycle source bank do i = 0, n_procs - 1 if (rank == i) then - ! ID's of first and last source particles - bank_first = rank*maxwork + 1 - bank_last = min((rank+1)*maxwork, n_particles) - - ! number of particles for this processor - work = bank_last - bank_first + 1 - do j = 1, work id = bank_first + j - 1 source_bank(j) % id = id @@ -216,4 +223,35 @@ contains end subroutine write_source_binary +!=============================================================================== +! READ_SOURCE reads a source distribution from a source.binary file and +! initializes the source bank +!=============================================================================== + + subroutine read_source_binary() + + integer(8) :: i ! index in source_bank + integer(8) :: n_sites ! number of sites in binary file + + ! Open binary source file for reading + open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='old', & + ACCESS='stream') + + ! Read number of source sites in file + read(UNIT=UNIT_SOURCE) n_sites + + do i = 1, maxwork + ! Set ID for source site + source_bank(i) % id = i ! bank_first + i - 1 + + ! Read position, angle, and energy + read(UNIT=UNIT_SOURCE) source_bank(i) % xyz, & + source_bank(i) % uvw, source_bank(i) % E + end do + + ! Close binary source file + close(UNIT=UNIT_SOURCE) + + end subroutine read_source_binary + end module source From 39603156ffd79bb76879e46cc6acce23edff4209 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 17:18:45 -0400 Subject: [PATCH 29/42] Can read binary source in parallel now. --- src/bank_header.F90 | 1 + src/source.F90 | 46 +++++++++++++++++++++++++++++++++++++++- src/utils/read_source.py | 2 +- 3 files changed, 47 insertions(+), 2 deletions(-) diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 0b42aa155a..dd255c5f0b 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -9,6 +9,7 @@ module bank_header !=============================================================================== type Bank + sequence integer(8) :: id ! Unique ID real(8) :: xyz(3) ! location of bank particle real(8) :: uvw(3) ! diretional cosines diff --git a/src/source.F90 b/src/source.F90 index a9de74bed1..f72bb76c54 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -11,6 +11,10 @@ module source use random_lcg, only: prn, set_particle_seed use string, only: to_str +#ifdef MPI + use mpi +#endif + implicit none contains @@ -232,6 +236,45 @@ contains integer(8) :: i ! index in source_bank integer(8) :: n_sites ! number of sites in binary file +#ifdef MPI + integer :: fh ! file handle + integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) +#endif + +#ifdef MPI + ! ========================================================================== + ! PARALLEL I/O USING MPI-2 ROUTINES + + ! Open binary source file for reading + call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_RDONLY, & + MPI_INFO_NULL, fh, mpi_err) + + offset = 0 + call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpi_err) + + ! Set proper offset for source data on this processor + offset = 8*(1 + rank*maxwork*7) + + do i = 1, work + ! Read position, angle, and energy -- note that we have used some + ! trickery here! The third argument gives the starting memory location, + ! and we have told it to read seven real(8)'s. Since xyz is only three + ! real(8)'s, it will write to the uvw and E members since those come next + ! in memory. We have guaranteed this by making the Bank type sequential. + + call MPI_FILE_READ_AT(fh, offset, source_bank(i) % xyz, 7, MPI_REAL8, & + MPI_STATUS_IGNORE, mpi_err) + offset = offset + 56 + + end do + + ! Close binary source file + call MPI_FILE_CLOSE(fh, mpi_err) + +#else + ! ========================================================================== + ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES ! Open binary source file for reading open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='old', & @@ -240,7 +283,7 @@ contains ! Read number of source sites in file read(UNIT=UNIT_SOURCE) n_sites - do i = 1, maxwork + do i = 1, work ! Set ID for source site source_bank(i) % id = i ! bank_first + i - 1 @@ -251,6 +294,7 @@ contains ! Close binary source file close(UNIT=UNIT_SOURCE) +#endif end subroutine read_source_binary diff --git a/src/utils/read_source.py b/src/utils/read_source.py index 6da8ff5008..dc43be9dfa 100755 --- a/src/utils/read_source.py +++ b/src/utils/read_source.py @@ -8,7 +8,7 @@ class SourceFile(object): self.f = open(filename, 'r') # Read number of source sites - self.n_sites = struct.unpack('l', self.f.read(8))[0] + self.n_sites = struct.unpack('q', self.f.read(8))[0] # Create list to store source sites self.sites = [] From 8dd7a57c33bf0610c313d2394c3f01af50243fa5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 17:54:52 -0400 Subject: [PATCH 30/42] Modified writing/reading of source binary (include 'id' attribute to take advantage of MPI_BANK). --- src/source.F90 | 37 ++++++++----------------------------- src/utils/read_source.py | 1 + 2 files changed, 9 insertions(+), 29 deletions(-) diff --git a/src/source.F90 b/src/source.F90 index f72bb76c54..ab5719e20c 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -207,8 +207,7 @@ contains subroutine write_source_binary() - integer(8) :: i - + ! open binary source file for writing open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', & ACCESS='stream') @@ -216,11 +215,7 @@ contains write(UNIT=UNIT_SOURCE) n_particles ! write information from the source bank - do i = 1_8, work - write(UNIT=UNIT_SOURCE) source_bank(i) % xyz - write(UNIT=UNIT_SOURCE) source_bank(i) % uvw - write(UNIT=UNIT_SOURCE) source_bank(i) % E - end do + write(UNIT=UNIT_SOURCE) source_bank(1:work) ! close binary source file close(UNIT=UNIT_SOURCE) @@ -234,7 +229,6 @@ contains subroutine read_source_binary() - integer(8) :: i ! index in source_bank integer(8) :: n_sites ! number of sites in binary file #ifdef MPI integer :: fh ! file handle @@ -254,20 +248,11 @@ contains MPI_STATUS_IGNORE, mpi_err) ! Set proper offset for source data on this processor - offset = 8*(1 + rank*maxwork*7) + offset = 8*(1 + rank*maxwork*8) - do i = 1, work - ! Read position, angle, and energy -- note that we have used some - ! trickery here! The third argument gives the starting memory location, - ! and we have told it to read seven real(8)'s. Since xyz is only three - ! real(8)'s, it will write to the uvw and E members since those come next - ! in memory. We have guaranteed this by making the Bank type sequential. - - call MPI_FILE_READ_AT(fh, offset, source_bank(i) % xyz, 7, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - offset = offset + 56 - - end do + ! Read all source sites + call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpi_err) ! Close binary source file call MPI_FILE_CLOSE(fh, mpi_err) @@ -283,14 +268,8 @@ contains ! Read number of source sites in file read(UNIT=UNIT_SOURCE) n_sites - do i = 1, work - ! Set ID for source site - source_bank(i) % id = i ! bank_first + i - 1 - - ! Read position, angle, and energy - read(UNIT=UNIT_SOURCE) source_bank(i) % xyz, & - source_bank(i) % uvw, source_bank(i) % E - end do + ! Read position, angle, and energy + read(UNIT=UNIT_SOURCE) source_bank(1:work) ! Close binary source file close(UNIT=UNIT_SOURCE) diff --git a/src/utils/read_source.py b/src/utils/read_source.py index dc43be9dfa..a7cb801c7f 100755 --- a/src/utils/read_source.py +++ b/src/utils/read_source.py @@ -15,6 +15,7 @@ class SourceFile(object): for i in range(self.n_sites): # Read position, angle, and energy + uid = struct.unpack('q', self.f.read(8)) xyz = struct.unpack('3d', self.f.read(24)) uvw = struct.unpack('3d', self.f.read(24)) E = struct.unpack('d', self.f.read(8)) From f1d9cef124e256cfbe2971159d75c704add8bbad Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 18:02:50 -0400 Subject: [PATCH 31/42] Reading/writing binary source now works in parallel. --- src/source.F90 | 40 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/src/source.F90 b/src/source.F90 index ab5719e20c..fb35728f7e 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -207,18 +207,50 @@ contains subroutine write_source_binary() - ! open binary source file for writing +#ifdef MPI + integer :: fh ! file handle + integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) + + ! ========================================================================== + ! PARALLEL I/O USING MPI-2 ROUTINES + + ! Open binary source file for reading + call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_CREATE + & + MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err) + + if (master) then + offset = 0 + call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpi_err) + end if + + ! Set proper offset for source data on this processor + offset = 8*(1 + rank*maxwork*8) + + ! Write all source sites + call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpi_err) + + ! Close binary source file + call MPI_FILE_CLOSE(fh, mpi_err) + +#else + ! ========================================================================== + ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES + + ! Open binary source file for writing open(UNIT=UNIT_SOURCE, FILE='source.binary', STATUS='replace', & ACCESS='stream') - ! write the number of particles + ! Write the number of particles write(UNIT=UNIT_SOURCE) n_particles - ! write information from the source bank + ! Write information from the source bank write(UNIT=UNIT_SOURCE) source_bank(1:work) - ! close binary source file + ! Close binary source file close(UNIT=UNIT_SOURCE) +#endif end subroutine write_source_binary From b940ec69f03e40c3738f42db03397b91ba4a5599 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 18:16:39 -0400 Subject: [PATCH 32/42] Fixed error in reading tag without coefficients. --- src/input_xml.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index d9385fd376..95199f319f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -69,6 +69,7 @@ contains energy_grid_ = "union" seed_ = 0_8 write_source = "" + source_ % type = "" ! Parse settings.xml file call read_xml_file_settings_t(filename) @@ -133,7 +134,7 @@ contains ! Verbosity if (verbosity_ > 0) verbosity = verbosity_ - if (associated(source_ % coeffs)) then + if (len(source_ % type) > 0) then ! Determine external source type type = source_ % type call lower_case(type) From cfc1dc2ab33f28c54013847373bd0eb4c17350f1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 18:20:12 -0400 Subject: [PATCH 33/42] Added check for source file size different from n_particles. --- src/source.F90 | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/source.F90 b/src/source.F90 index fb35728f7e..45e9d76996 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -275,10 +275,18 @@ contains call MPI_FILE_OPEN(MPI_COMM_WORLD, 'source.binary', MPI_MODE_RDONLY, & MPI_INFO_NULL, fh, mpi_err) + ! Read number of source sites in file offset = 0 call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, & MPI_STATUS_IGNORE, mpi_err) + ! Check that number of source sites matches + if (n_sites /= n_particles) then + message = "No support yet for source files of different size than & + &specified number of particles per generation." + call fatal_error() + end if + ! Set proper offset for source data on this processor offset = 8*(1 + rank*maxwork*8) @@ -300,6 +308,13 @@ contains ! Read number of source sites in file read(UNIT=UNIT_SOURCE) n_sites + ! Check that number of source sites matches + if (n_sites /= n_particles) then + message = "No support yet for source files of different size than & + &specified number of particles per generation." + call fatal_error() + end if + ! Read position, angle, and energy read(UNIT=UNIT_SOURCE) source_bank(1:work) From aee6beb01a48d237ecfe2e66774f9cfdd60c89cf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 22 Mar 2012 18:28:40 -0400 Subject: [PATCH 34/42] Fixed bug in input_xml causing gfortran to not compile. --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 95199f319f..e334007fda 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -68,7 +68,7 @@ contains verbosity_ = 0 energy_grid_ = "union" seed_ = 0_8 - write_source = "" + write_source_ = "" source_ % type = "" ! Parse settings.xml file From 2c24c3d4057eb9d53e56939fc706655a553f5e8f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 23 Mar 2012 08:42:03 -0400 Subject: [PATCH 35/42] Added git SHA1 hash on output through Makefile symbol. --- src/Makefile | 8 +++++++- src/output.F90 | 3 +++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/Makefile b/src/Makefile index 14b74bd85a..a5c9237e2f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -22,6 +22,12 @@ OPTIMIZE = no USE_MPI = no USE_HDF5 = no +#=============================================================================== +# Add git SHA-1 hash +#=============================================================================== + +GIT_SHA1 = $(shell git log -1 | head -n 1 | awk '{print $$2}') + #=============================================================================== # GNU Fortran compiler options #=============================================================================== @@ -212,7 +218,7 @@ neat: .PHONY: all xml-fortran clean neat distclean %.o: %.F90 - $(F90) $(F90FLAGS) -Ixml-fortran -Ixml-fortran/templates -c $< + $(F90) $(F90FLAGS) -DGIT_SHA1="\"$(GIT_SHA1)\"" -Ixml-fortran -Ixml-fortran/templates -c $< #=============================================================================== # Dependencies diff --git a/src/output.F90 b/src/output.F90 index 920d9643f7..c389fab8d1 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -50,6 +50,9 @@ contains ' Developed At: Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +#ifdef GIT_SHA1 + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 +#endif ! Write the date and time call get_today(today_date, today_time) From 0adb96d042078548bb01de8b915aa9bd919bee0e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Mar 2012 11:19:39 -0400 Subject: [PATCH 36/42] Cleaned up source, added comments, removed unnecessary variables. --- src/DEPENDENCIES | 4 +- src/ace.F90 | 7 +- src/bank_header.F90 | 4 + src/energy_grid.F90 | 2 +- src/initialize.F90 | 14 ++-- src/input_xml.F90 | 3 +- src/main.F90 | 1 - src/mesh_header.F90 | 14 ++-- src/output.F90 | 190 ++++++++++++++++++++++-------------------- src/physics.F90 | 76 +---------------- src/plot.F90 | 111 +++++++++++++----------- src/plot_header.F90 | 2 +- src/ppmlib.F90 | 66 ++++++++++----- src/search.F90 | 1 - src/source.F90 | 2 +- src/source_header.F90 | 2 +- src/tally.F90 | 6 +- src/timing.F90 | 10 +-- 18 files changed, 245 insertions(+), 270 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index bcc3b533ff..d43e0f3539 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -99,6 +99,7 @@ initialize.o: random_lcg.o initialize.o: source.o initialize.o: string.o initialize.o: tally.o +initialize.o: tally_header.o initialize.o: timing.o input_xml.o: constants.o @@ -143,7 +144,6 @@ main.o: intercycle.o main.o: output.o main.o: physics.o main.o: plot.o -main.o: random_lcg.o main.o: source.o main.o: string.o main.o: tally.o @@ -177,6 +177,7 @@ physics.o: geometry.o physics.o: geometry_header.o physics.o: global.o physics.o: interpolation.o +physics.o: material_header.o physics.o: output.o physics.o: particle_header.o physics.o: random_lcg.o @@ -198,7 +199,6 @@ plot.o: ppmlib.o plot.o: source.o plot.o: string.o -search.o: constants.o search.o: error.o search.o: global.o diff --git a/src/ace.F90 b/src/ace.F90 index 0406c7dbca..af256c9e21 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -22,8 +22,6 @@ module ace real(8), allocatable :: XSS(:) ! Cross section data integer :: XSS_index ! current index in XSS data - type(DictionaryCI), pointer :: already_read => null() - private :: NXS private :: JXS private :: XSS @@ -47,6 +45,7 @@ contains type(Material), pointer :: mat => null() type(Nuclide), pointer :: nuc => null() type(SAB_Table), pointer :: sab => null() + type(DictionaryCI), pointer :: already_read => null() ! ========================================================================== ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES @@ -154,6 +153,8 @@ contains name = xs_listings(index_list) % name alias = xs_listings(index_list) % alias + ! Read the ACE table into the appropriate entry on the nuclides + ! array call read_ace_table(index_nuclides, index_list) ! Print out information on table to cross_sections.out file @@ -172,6 +173,8 @@ contains index_list = dict_get_key(xs_listing_dict, name) index_sab = dict_get_key(sab_dict, name) + ! Read the ACE table into the appropriate entry on the sab_tables + ! array call read_ace_table(index_sab, index_list) call dict_add_key(already_read, name, 0) diff --git a/src/bank_header.F90 b/src/bank_header.F90 index dd255c5f0b..233ff611d4 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -9,7 +9,11 @@ module bank_header !=============================================================================== type Bank + ! The 'sequence' attribute is used here to ensure that the data listed + ! appears in the given order. This is important for MPI purposes when bank + ! sites are sent from one processor to another. sequence + integer(8) :: id ! Unique ID real(8) :: xyz(3) ! location of bank particle real(8) :: uvw(3) ! diretional cosines diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 3fb1d940d2..4828338a2e 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -37,7 +37,7 @@ contains ! Create dictionary for keeping track of cross sections already added call dict_create(already_added) - ! Loop over all files + ! Loop over all materials do i = 1, n_materials mat => materials(i) diff --git a/src/initialize.F90 b/src/initialize.F90 index 6bc2e899c8..c119cd9bab 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -3,13 +3,13 @@ module initialize use ace, only: read_xs use bank_header, only: Bank use constants - use datatypes, only: dict_create, dict_add_key, dict_get_key, & - dict_has_key, dict_keys + use datatypes, only: dict_create, dict_get_key, dict_has_key, & + dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII use energy_grid, only: unionized_grid use error, only: fatal_error use geometry, only: neighbor_lists - use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE + use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml @@ -19,7 +19,8 @@ module initialize use random_lcg, only: initialize_prng use source, only: initialize_source use string, only: to_str, starts_with, ends_with, lower_case - use tally, only: create_tally_map, TallyObject + use tally, only: create_tally_map + use tally_header, only: TallyObject use timing, only: timer_start, timer_stop #ifdef MPI @@ -33,8 +34,6 @@ module initialize implicit none - type(DictionaryII), pointer :: build_dict => null() - contains !=============================================================================== @@ -286,9 +285,6 @@ contains ! Create special dictionary used in input_xml call dict_create(cells_in_univ_dict) - - ! Create special dictionary for building universes - call dict_create(build_dict) end subroutine create_dictionaries diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e334007fda..039f8cb653 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1,8 +1,7 @@ module input_xml use constants - use datatypes, only: dict_create, dict_add_key, dict_has_key, & - dict_get_key + use datatypes, only: dict_add_key, dict_has_key, dict_get_key use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice use global diff --git a/src/main.F90 b/src/main.F90 index 31092f79f6..bf5b922bc4 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -8,7 +8,6 @@ program main use output, only: write_message, header use plotter, only: run_plot use physics, only: transport - use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str use tally, only: synchronize_tallies, add_to_score diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index b70bbd45b4..8088972ce0 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -8,13 +8,13 @@ module mesh_header !=============================================================================== type StructuredMesh - integer :: id - integer :: type - integer :: n_dimension - integer, allocatable :: dimension(:) - real(8), allocatable :: lower_left(:) - real(8), allocatable :: upper_right(:) - real(8), allocatable :: width(:) + integer :: id ! user-specified id + integer :: type ! rectangular, hexagonal + integer :: n_dimension ! rank of mesh + integer, allocatable :: dimension(:) ! number of cells in each direction + real(8), allocatable :: lower_left(:) ! lower-left corner of mesh + real(8), allocatable :: upper_right(:) ! upper-right corner of mesh + real(8), allocatable :: width(:) ! width of each mesh cell end type StructuredMesh end module mesh_header diff --git a/src/output.F90 b/src/output.F90 index c389fab8d1..2fc6c77e96 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -29,8 +29,8 @@ contains subroutine title() - character(10) :: today_date - character(8) :: today_time + character(10) :: date_ + character(10) :: time_ write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & ' .d88888b. 888b d888 .d8888b.', & @@ -55,18 +55,23 @@ contains #endif ! Write the date and time - call get_today(today_date, today_time) + call date_and_time(DATE=date_, TIME=time_) + date_ = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) + time_ = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(time_) ! Write information to summary file call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011 Massachusetts Institute of Technology" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +#ifdef GIT_SHA1 + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 +#endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(time_) ! Write information on number of processors #ifdef MPI @@ -89,8 +94,8 @@ contains integer, optional :: unit ! unit to write to integer, optional :: level ! specified header level - integer :: n - integer :: m + integer :: n ! number of = signs on left + integer :: m ! number of = signs on right integer :: unit_ ! unit to write to integer :: header_level ! actual header level character(MAX_LINE_LEN) :: line @@ -140,14 +145,15 @@ contains subroutine write_message(level) - integer, optional :: level + integer, optional :: level ! verbosity level - integer :: n_lines - integer :: i + integer :: n_lines ! number of lines needed + integer :: i ! index for lines ! Only allow master to print to screen if (.not. master .and. present(level)) return + ! TODO: Take care of line wrapping so words don't get cut off if (.not. present(level) .or. level <= verbosity) then n_lines = (len_trim(message)-1)/79 + 1 do i = 1, n_lines @@ -157,55 +163,13 @@ contains end subroutine write_message -!=============================================================================== -! GET_TODAY determines the date and time at which the program began execution -! and returns it in a readable format -!=============================================================================== - - subroutine get_today(today_date, today_time) - - character(10), intent(out) :: today_date - character(8), intent(out) :: today_time - - integer :: val(8) - character(8) :: date_ - character(10) :: time_ - character(5) :: zone - - call date_and_time(date_, time_, zone, val) - ! val(1) = year (YYYY) - ! val(2) = month (MM) - ! val(3) = day (DD) - ! val(4) = timezone - ! val(5) = hours (HH) - ! val(6) = minutes (MM) - ! val(7) = seconds (SS) - ! val(8) = milliseconds - - if (val(2) < 10) then - if (val(3) < 10) then - today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(6:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - else - if (val(3) < 10) then - today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(5:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - end if - today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) - - end subroutine get_today - !=============================================================================== ! PRINT_PARTICLE displays the attributes of a particle !=============================================================================== subroutine print_particle() - integer :: i + integer :: i ! index for coordinate levels type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() @@ -266,6 +230,7 @@ contains write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface)) end if + ! Display weight, energy, grid index, and interpolation factor write(ou,*) ' Weight = ' // to_str(p % wgt) write(ou,*) ' Energy = ' // to_str(p % E) write(ou,*) ' IE = ' // to_str(p % IE) @@ -301,27 +266,35 @@ contains subroutine print_cell(c, unit) type(Cell), pointer :: c - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: temp - integer :: i - integer :: unit_ + integer :: index_cell ! index in cells array + integer :: i ! loop index for surfaces + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() + ! Set unit to stdout if not already set if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id for cell write(unit_,*) 'Cell ' // to_str(c % id) - temp = dict_get_key(cell_dict, c % id) - write(unit_,*) ' Array Index = ' // to_str(temp) + + ! Find index in cells array and write + index_cell = dict_get_key(cell_dict, c % id) + write(unit_,*) ' Array Index = ' // to_str(index_cell) + + ! Write what universe this cell is in u => universes(c % universe) write(unit_,*) ' Universe = ' // to_str(u % id) + + ! Write information on fill for cell select case (c % type) case (CELL_NORMAL) write(unit_,*) ' Fill = NONE' @@ -332,12 +305,16 @@ contains l => lattices(c % fill) write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select + + ! Write information on material if (c % material == 0) then write(unit_,*) ' Material = NONE' else m => materials(c % material) write(unit_,*) ' Material = ' // to_str(m % id) end if + + ! Write surface specification string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -367,24 +344,31 @@ contains type(Universe), pointer :: univ integer, optional :: unit - integer :: i - integer :: unit_ + integer :: i ! loop index for cells in this universe + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() - type(Universe), pointer :: base_u + type(Universe), pointer :: base_u => null() + ! Set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Get a pointer to the base universe base_u => universes(BASE_UNIVERSE) + ! Write user-specified id for this universe write(unit_,*) 'Universe ' // to_str(univ % id) + + ! If this is the base universe, indicate so if (associated(univ, base_u)) then write(unit_,*) ' Base Universe' end if + + ! Write list of cells in this universe string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) @@ -404,14 +388,16 @@ contains type(Lattice), pointer :: lat integer, optional :: unit - integer :: unit_ + integer :: unit_ ! unit to write to + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write information about lattice write(unit_,*) 'Lattice ' // to_str(lat % id) write(unit_,*) ' n_x = ' // to_str(lat % n_x) write(unit_,*) ' n_y = ' // to_str(lat % n_y) @@ -430,19 +416,23 @@ contains subroutine print_surface(surf, unit) type(Surface), pointer :: surf - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: i - integer :: unit_ + integer :: i ! loop index for coefficients + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of surface write(unit_,*) 'Surface ' // to_str(surf % id) + + ! Write type of surface select case (surf % type) case (SURF_PX) string = "X Plane" @@ -469,12 +459,14 @@ contains end select write(unit_,*) ' Type = ' // trim(string) + ! Write coefficients for this surface string = "" do i = 1, size(surf % coeffs) string = trim(string) // ' ' // to_str(surf % coeffs(i), 4) end do write(unit_,*) ' Coefficients = ' // trim(string) + ! Write neighboring cells on positive side of this surface string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) @@ -483,6 +475,7 @@ contains end if write(unit_,*) ' Positive Neighbors = ' // trim(string) + ! Write neighboring cells on negative side of this surface string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) @@ -490,6 +483,8 @@ contains end do end if write(unit_,*) ' Negative Neighbors =' // trim(string) + + ! Write boundary condition for this surface select case (surf % bc) case (BC_TRANSMIT) write(unit_,*) ' Boundary Condition = Transmission' @@ -513,12 +508,13 @@ contains type(Material), pointer :: mat integer, optional :: unit - integer :: i - integer :: unit_ - real(8) :: density + integer :: i ! loop index for nuclides + integer :: unit_ ! unit to write to + real(8) :: density ! density in atom/b-cm character(MAX_LINE_LEN) :: string type(Nuclide), pointer :: nuc => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else @@ -559,9 +555,9 @@ contains type(TallyObject), pointer :: t integer, optional :: unit - integer :: i - integer :: id - integer :: unit_ + integer :: i ! index for filter or score bins + integer :: id ! user-specified id + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() @@ -569,14 +565,17 @@ contains type(Material), pointer :: m => null() type(StructuredMesh), pointer :: sm => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of tally write(unit_,*) 'Tally ' // to_str(t % id) + ! Write any cells bins if present if (t % n_filter_bins(FILTER_CELL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELL) @@ -587,6 +586,7 @@ contains write(unit_, *) ' Cell Bins:' // trim(string) end if + ! Write any surface bins if present if (t % n_filter_bins(FILTER_SURFACE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_SURFACE) @@ -597,6 +597,7 @@ contains write(unit_, *) ' Surface Bins:' // trim(string) end if + ! Write any universe bins if present if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_UNIVERSE) @@ -607,6 +608,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any material bins if present if (t % n_filter_bins(FILTER_MATERIAL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_MATERIAL) @@ -617,6 +619,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any mesh bins if present if (t % n_filter_bins(FILTER_MESH) > 0) then string = "" id = t % mesh @@ -628,6 +631,7 @@ contains write(unit_, *) ' Mesh Bins:' // trim(string) end if + ! Write any birth region bins if present if (t % n_filter_bins(FILTER_CELLBORN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELLBORN) @@ -638,6 +642,7 @@ contains write(unit_, *) ' Birth Region Bins:' // trim(string) end if + ! Write any incoming energy bins if present if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1 @@ -647,6 +652,7 @@ contains write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if + ! Write any outgoing energy bins if present if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1 @@ -656,6 +662,7 @@ contains write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if + ! Write any score bins if present if (t % n_score_bins > 0) then string = "" do i = 1, t % n_score_bins @@ -674,7 +681,7 @@ contains string = trim(string) // ' nu-fission' end select end do - write(unit_,*) ' Macro Reactions:' // trim(string) + write(unit_,*) ' Scores:' // trim(string) end if write(unit_,*) @@ -687,7 +694,7 @@ contains subroutine print_geometry() - integer :: i + integer :: i ! loop index for various arrays type(Surface), pointer :: s => null() type(Cell), pointer :: c => null() type(Universe), pointer :: u => null() @@ -735,12 +742,12 @@ contains type(Nuclide), pointer :: nuc integer, optional :: unit - integer :: i - integer :: unit_ - integer :: size_total - integer :: size_xs - integer :: size_angle - integer :: size_energy + integer :: i ! loop index over nuclides + integer :: unit_ ! unit to write to + integer :: size_total ! memory used by nuclide (bytes) + integer :: size_xs ! memory used for cross-sections (bytes) + integer :: size_angle ! memory used for angle distributions (bytes) + integer :: size_energy ! memory used for energy distributions (bytes) type(Reaction), pointer :: rxn => null() type(UrrData), pointer :: urr => null() @@ -824,7 +831,7 @@ contains subroutine print_summary() - integer :: i + integer :: i ! loop index character(15) :: string type(Material), pointer :: m => null() type(TallyObject), pointer :: t => null() @@ -896,24 +903,24 @@ contains subroutine print_plot() - integer i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() ! Display header for plotting call header("PLOTTING SUMMARY") - do i=1,n_plots + do i = 1, n_plots pl => plots(i) - ! Print plot id + ! Write plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) - ! Print plotting origin + ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & " " // trim(to_str(pl % origin(2))) // " " // & trim(to_str(pl % origin(3))) - ! Print plotting width + ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & @@ -940,8 +947,8 @@ contains subroutine print_runtime() - integer(8) :: total_particles - real(8) :: speed + integer(8) :: total_particles ! total # of particles simulated + real(8) :: speed ! # of neutrons/second character(15) :: string ! display header block @@ -1017,7 +1024,8 @@ contains end subroutine create_summary_file !=============================================================================== -! CREATE_XS_SUMMARY_FILE +! CREATE_XS_SUMMARY_FILE creates an output file to write information about the +! cross section tables used in the simulation. !=============================================================================== subroutine create_xs_summary_file() diff --git a/src/physics.F90 b/src/physics.F90 index 27ec0adc20..a36202ea72 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -3,7 +3,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants use cross_section, only: calculate_xs - use endf, only: reaction_name, is_fission, is_scatter + use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_prompt, nu_delayed use geometry, only: find_cell, distance_to_boundary, cross_surface, & @@ -11,6 +11,7 @@ module physics use geometry_header, only: Universe, BASE_UNIVERSE use global use interpolation, only: interpolate_tab1 + use material_header, only: Material use output, only: write_message use particle_header, only: LocalCoord use random_lcg, only: prn @@ -2061,77 +2062,4 @@ contains end function watt_spectrum -!=============================================================================== -! WIGNER samples a Wigner distribution of energy level spacings. Note that this -! scheme is C50 in the Monte Carlo Sampler from Los Alamos (LA-9721-MS). -!=============================================================================== - - function wigner(D_avg) result (D) - - real(8), intent(in) :: D_avg ! average level spacing - real(8) :: D ! sampled level spacing - - real(8) :: c - - c = -4.*D_avg*D_avg/PI * log(prn()) - D = sqrt(c) - - end function wigner - -!=============================================================================== -! CHI_SQUARED samples a chi-squared distribution with n degrees of freedom. The -! distribution of resonance widths in the unresolved region is given by a -! chi-squared distribution. For the special case of n=1, this is a Porter-Thomas -! distribution. For cases with n odd, rule C64 is used whereas for cases with n -! even, rule C45 is used. -!=============================================================================== - - function chi_squared(n, G_avg) result(G) - - integer, intent(in) :: n ! number of degrees of freedom - real(8), intent(in), optional :: G_avg ! average resonance width - - integer :: i ! loop index - real(8) :: G ! sampled random variable (or resonance width) - real(8) :: x, y, c ! dummy variables - real(8) :: r1, r2 ! psuedorandom numbers - - select case (mod(n,2)) - case (0) - ! Even number of degrees of freedom can be sampled via rule C45. We can - ! sample x as -2/n*log(product(r_i, i = 1 to n/2)) - x = ONE - do i = 1, n/2 - x = x * prn() - end do - x = -2./n * log(x) - - case (1) - ! Odd number of degrees of freedom can be sampled via rule C64. We can - ! sample x as -2/n*(log(r)*cos^2(pi/2*r) + log(product(r_i, i = 1 to - ! floor(n/2))) - - ! Note that we take advantage of integer division on n/2 - y = ONE - do i = 1, n/2 - y = y * prn() - end do - - r1 = prn() - r2 = prn() - c = cos(PI/2.*r2) - x = -2./n * (log(y) + log(r1)*c*c) - end select - - ! If sampling a chi-squared distribution for a resonance width and the - ! average resonance width has been given, return the sampled resonance - ! width. - if (present(G_avg)) then - G = x * G_avg - else - G = x - end if - - end function chi_squared - end module physics diff --git a/src/plot.F90 b/src/plot.F90 index 61d0ea74dd..7582130454 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -2,12 +2,11 @@ module plotter use constants use error, only: fatal_error - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice, cell_contains + use geometry, only: find_cell use geometry_header, only: Universe, BASE_UNIVERSE use global use output, only: write_message - use particle_header, only: LocalCoord, deallocate_coord + use particle_header, only: deallocate_coord use plot_header use ppmlib, only: Image, init_image, allocate_image, & deallocate_image, set_pixel @@ -18,17 +17,16 @@ module plotter contains - !=============================================================================== -! RUN_PLOT +! RUN_PLOT controls the logic for making one or many plots !=============================================================================== subroutine run_plot() - integer :: i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() - do i=1,n_plots + do i = 1, n_plots pl => plots(i) ! Display output message @@ -36,31 +34,34 @@ contains call write_message(5) if (pl % type == PLOT_TYPE_SLICE) then - ! create 2d image call create_ppm(pl) - end if end do end subroutine run_plot !=============================================================================== -! create_ppm +! CREATE_PPM creates an image based on user input from a plots.xml +! specification in the portable pixmap format (PPM) !=============================================================================== + subroutine create_ppm(pl) - type(Plot), pointer :: pl + type(Plot), pointer :: pl - type(Image) :: img - integer :: in_i, out_i - integer :: x, y - integer :: r, g, b - real(8) :: in_pixel, out_pixel - real(8) :: xyz(3) - logical :: found_cell - type(Cell), pointer :: c => null() + integer :: in_i + integer :: out_i + integer :: x, y ! pixel location + integer :: r, g, b ! colors (red, green, blue) from 0-255 + real(8) :: in_pixel + real(8) :: out_pixel + real(8) :: xyz(3) + logical :: found_cell + type(Image) :: img + type(Cell), pointer :: c => null() + ! Initialize and allocate space for image call init_image(img) call allocate_image(img, pl % pixels(1), pl % pixels(2)) @@ -97,8 +98,8 @@ contains p % coord % uvw = (/ 1, 0, 0 /) p % coord % universe = BASE_UNIVERSE - do y=1, img%height - do x=1, img%width + do y = 1, img % height + do x = 1, img % width call deallocate_coord(p % coord0 % next) p % coord => p % coord0 @@ -106,63 +107,74 @@ contains call find_cell(found_cell) if (.not. found_cell) then - r = pl % not_found % rgb(1) - g = pl % not_found % rgb(2) - b = pl % not_found % rgb(3) + ! If no cell, revert to default color + r = pl % not_found % rgb(1) + g = pl % not_found % rgb(2) + b = pl % not_found % rgb(3) else - if (pl % color_by == PLOT_COLOR_MATS) then - c => cells(p % coord % cell) - r = pl % colors(c % material) % rgb(1) - g = pl % colors(c % material) % rgb(2) - b = pl % colors(c % material) % rgb(3) - else if (pl % color_by == PLOT_COLOR_CELLS) then - r = pl % colors(p % coord % cell) % rgb(1) - g = pl % colors(p % coord % cell) % rgb(2) - b = pl % colors(p % coord % cell) % rgb(3) - else - r = 0 - g = 0 - b = 0 - end if + if (pl % color_by == PLOT_COLOR_MATS) then + ! Assign color based on material + c => cells(p % coord % cell) + r = pl % colors(c % material) % rgb(1) + g = pl % colors(c % material) % rgb(2) + b = pl % colors(c % material) % rgb(3) + else if (pl % color_by == PLOT_COLOR_CELLS) then + ! Assign color based on cell + r = pl % colors(p % coord % cell) % rgb(1) + g = pl % colors(p % coord % cell) % rgb(2) + b = pl % colors(p % coord % cell) % rgb(3) + else + r = 0 + g = 0 + b = 0 + end if end if + ! Create a pixel at (x,y) with color (r,g,b) call set_pixel(img, x, y, r, g, b) + ! Advance pixel in first direction p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel end do + ! Advance pixel in second direction p % coord0 % xyz(in_i) = xyz(in_i) p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel end do + ! Write out the ppm to a file call output_ppm(pl,img) + ! Free up space call deallocate_image(img) end subroutine create_ppm - !=============================================================================== -! output_ppm +! OUTPUT_PPM writes out a previously generated image to a PPM file !=============================================================================== - subroutine output_ppm(pl,img) - type(Plot), pointer :: pl - type(Image), intent(in) :: img + subroutine output_ppm(pl, img) - integer :: i, j + type(Plot), pointer :: pl + type(Image), intent(in) :: img + integer :: i ! loop index for height + integer :: j ! loop index for width + + ! Open PPM file for writing open(UNIT=UNIT_PLOT, FILE=pl % path_plot) + ! Write header write(UNIT_PLOT, '(A2)') 'P6' write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height write(UNIT_PLOT, '(A)') '255' - do j=1, img%height - do i=1, img%width + ! 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)) + achar(img%green(i,j)), achar(img%blue(i,j)) end do end do @@ -171,5 +183,4 @@ contains end subroutine output_ppm - end module plotter diff --git a/src/plot_header.F90 b/src/plot_header.F90 index c0ae470643..6268ffadf9 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -26,7 +26,7 @@ module plot_header integer :: basis ! direction of plot slice integer :: pixels(2) ! 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(ObjectColor), allocatable :: colors(:) ! colors of cells/mats end type Plot integer :: PLOT_TYPE_SLICE = 1 diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 index 2c9687617a..d67268412e 100644 --- a/src/ppmlib.F90 +++ b/src/ppmlib.F90 @@ -13,6 +13,10 @@ module ppmlib contains +!=============================================================================== +! INIT_IMAGE initializes the Image derived type +!=============================================================================== + subroutine init_image(img) type(Image), intent(out) :: img @@ -26,29 +30,46 @@ contains 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, h + 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) + 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) @@ -58,13 +79,15 @@ contains inside = .false. - if ( ( x < img % width) .and. & - ( y < img % height ) .and. & - ( x >= 0 ) .and. & - ( y >= 0 ) ) inside = .true. + 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) @@ -73,29 +96,32 @@ contains valid = .false. - if ( img % width == 0 ) return - if ( img % height == 0 ) return - if ( .not. associated(img % red) .or. & + 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 + .not. associated(img % blue)) return valid = .true. end function valid_image - +!=============================================================================== +! SET_PIXEL sets the +!=============================================================================== + subroutine set_pixel(img, x, y, r, g, b) type(Image), intent(inout) :: img - integer, intent(in) :: x, y - integer, intent(in) :: r, g, b + 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) + 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 subroutine set_pixel end module ppmlib diff --git a/src/search.F90 b/src/search.F90 index 3cdb51bfcb..5968113fec 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -1,6 +1,5 @@ module search - use constants, only: ONE, MAX_LINE_LEN use error, only: fatal_error use global, only: message diff --git a/src/source.F90 b/src/source.F90 index 45e9d76996..81d6ab2e6d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,7 +1,7 @@ module source use bank_header, only: Bank - use constants, only: ONE + use constants use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE use global diff --git a/src/source_header.F90 b/src/source_header.F90 index 077e6338d8..638a960244 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -8,7 +8,7 @@ module source_header !=============================================================================== type ExtSource - integer :: type ! type of source, e.g. 'box' or 'cell' + integer :: type ! type of source, e.g. 'box' or 'point' real(8), allocatable :: values(:) ! values for particular source type end type ExtSource diff --git a/src/tally.F90 b/src/tally.F90 index db57a56802..50d29e6061 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1062,7 +1062,8 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_BATCH_ESTIMATE +! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many +! generations) into a single realization of a random variable. !=============================================================================== elemental subroutine accumulate_batch_estimate(score) @@ -1540,7 +1541,8 @@ contains end function get_label !=============================================================================== -! CALCULATE_STATISTICS +! CALCULATE_STATISTICS determines the sample mean and the standard deviation of +! the mean for a TallyScore. !=============================================================================== elemental subroutine calculate_statistics(score) diff --git a/src/timing.F90 b/src/timing.F90 index f6e91c0583..483057125d 100644 --- a/src/timing.F90 +++ b/src/timing.F90 @@ -13,7 +13,7 @@ module timing type Timer logical :: running = .false. ! is timer running? integer :: start_counts = 0 ! counts when started - real(8) :: elapsed = 0. ! total time elapsed in seconds + real(8) :: elapsed = ZERO ! total time elapsed in seconds !!$ contains !!$ procedure :: start => timer_start !!$ procedure :: get_value => timer_get_value @@ -24,7 +24,7 @@ module timing contains !=============================================================================== -! TIMER_START +! TIMER_START starts running a timer and measures the current time !=============================================================================== subroutine timer_start(self) @@ -38,7 +38,7 @@ contains end subroutine timer_start !=============================================================================== -! TIMER_GET_VALUE +! TIMER_GET_VALUE returns the current value of the timer !=============================================================================== function timer_get_value(self) result(elapsed) @@ -61,7 +61,7 @@ contains end function timer_get_value !=============================================================================== -! TIMER_STOP +! TIMER_STOP stops the timer and sets the elapsed time !=============================================================================== subroutine timer_stop(self) @@ -78,7 +78,7 @@ contains end subroutine timer_stop !=============================================================================== -! TIMER_RESET +! TIMER_RESET resets a timer to have a zero value !=============================================================================== subroutine timer_reset(self) From 93bfdcccdd2e23cc9c3f9fdd86ec216d295fdc0f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 25 Mar 2012 17:18:26 -0400 Subject: [PATCH 37/42] Logic for tracklength tallies with mesh filter. Answers still unreliable. --- src/mesh.F90 | 87 ++++++++++++++++++ src/tally.F90 | 239 +++++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 325 insertions(+), 1 deletion(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 756525327c..3e489903b4 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -217,4 +217,91 @@ contains end subroutine count_fission_sites +!=============================================================================== +! MESH_INTERSECT +!=============================================================================== + + function mesh_intersects(m, xyz0, xyz1) result(intersects) + + type(StructuredMesh), pointer :: m + real(8), intent(in) :: xyz0(3) + real(8), intent(in) :: xyz1(3) + logical :: intersects + + real(8) :: x0, y0, z0 ! track start point + real(8) :: x1, y1, z1 ! track end point + real(8) :: xi, yi, zi ! track intersection point with mesh + + ! Copy coordinates of starting point + x0 = xyz0(1) + y0 = xyz0(2) + z0 = xyz0(3) + + ! Copy coordinates of ending point + x1 = xyz1(1) + y1 = xyz1(2) + z1 = xyz1(3) + + ! Check if line intersects bottom surface -- calculate the intersection + ! point (y,z) + xi = m % lower_left(1) + yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) + if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + intersects = .true. + return + end if + + ! Check if line intersects left surface -- calculate the intersection point + ! (x,z) + yi = m % lower_left(2) + xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) + if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + intersects = .true. + return + end if + + ! Check if line intersects front surface -- calculate the intersection point + ! (x,y) + zi = m % lower_left(3) + xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) + if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + intersects = .true. + return + end if + + ! Check if line intersects top surface -- calculate the intersection + ! point (y,z) + xi = m % upper_right(1) + yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) + if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + intersects = .true. + return + end if + + ! Check if line intersects right surface -- calculate the intersection point + ! (x,z) + yi = m % upper_right(2) + xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) + if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + intersects = .true. + return + end if + + ! Check if line intersects back surface -- calculate the intersection point + ! (x,y) + zi = m % upper_right(3) + xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) + if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + intersects = .true. + return + end if + + end function mesh_intersects + end module mesh diff --git a/src/tally.F90 b/src/tally.F90 index 50d29e6061..1f89af152e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3,7 +3,8 @@ module tally use constants use error, only: fatal_error use global - use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices + use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & + mesh_indices_to_bin, mesh_intersects use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search @@ -547,6 +548,14 @@ contains do i = 1, n_tracklength_tallies t => tallies(tracklength_tallies(i)) + ! Check if this tally has a mesh filter -- if so, we treat it separately + ! since multiple bins can be scored to with a single track + + if (t % n_filter_bins(FILTER_MESH) > 0) then + call score_tl_on_mesh(tracklength_tallies(i)) + cycle + end if + ! ======================================================================= ! DETERMINE SCORING BIN COMBINATION @@ -606,6 +615,234 @@ contains end subroutine score_tracklength_tally +!=============================================================================== +! SCORE_TL_ON_MESH calculate fluxes and reaction rates based on the track-length +! estimate of the flux specifically for tallies that have mesh filters. For +! these tallies, it is possible to score to multiple mesh cells for each track. +!=============================================================================== + + subroutine score_tl_on_mesh(index_tally) + + integer, intent(in) :: index_tally + + integer :: i ! loop index for filter/score bins + integer :: j ! loop index for direction + integer :: k ! loop index for mesh cell crossings + integer :: ijk0(3) ! indices of starting coordinates + integer :: ijk1(3) ! indices of ending coordinates + integer :: ijk_cross(3) ! indices of mesh cell crossed + integer :: n_cross ! number of surface crossings + integer :: bins(N_FILTER_TYPES) ! scoring bin combination + integer :: score_index ! single index for single bin + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + real(8) :: flux ! tracklength estimate of flux + real(8) :: score ! actual score (e.g., flux*xs) + real(8) :: uvw(3) ! cosine of angle of particle + real(8) :: xyz0(3) ! starting/intermediate coordinates + real(8) :: xyz1(3) ! ending coordinates of particle + real(8) :: xyz_cross(3) ! coordinates of next boundary + real(8) :: d(3) ! distance to each bounding surface + real(8) :: distance ! distance traveled in mesh cell + logical :: found_bin ! was a scoring bin found? + logical :: start_in_mesh ! starting coordinates inside mesh? + logical :: end_in_mesh ! ending coordinates inside mesh? + type(TallyObject), pointer :: t => null() + type(StructuredMesh), pointer :: m => null() + + found_bin = .true. + bins = 1 + t => tallies(index_tally) + + ! ========================================================================== + ! CHECK IF THIS TRACK INTERSECTS THE MESH + + ! Copy starting and ending location of particle + xyz0 = p % last_xyz + xyz1 = p % coord0 % xyz + + ! Determine indices for starting and ending location + m => meshes(t % mesh) + call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) + call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) + + ! Check if start or end is in mesh -- if not, check if track still + ! intersects with mesh + if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + if (.not. mesh_intersects(m, xyz0, xyz1)) return + end if + + ! ========================================================================= + ! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH + + FILTER_LOOP: do i = 1, t % n_filters + + select case (t % filters(i)) + case (FILTER_UNIVERSE) + ! determine next universe bin + ! TODO: Account for multiple universes when performing this filter + bins(FILTER_UNIVERSE) = get_next_bin(FILTER_UNIVERSE, & + p % coord % universe, index_tally) + if (bins(FILTER_UNIVERSE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_MATERIAL) + bins(FILTER_MATERIAL) = get_next_bin(FILTER_MATERIAL, & + p % material, index_tally) + if (bins(FILTER_MATERIAL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELL) + ! determine next cell bin + ! TODO: Account for cells in multiple levels when performing this filter + bins(FILTER_CELL) = get_next_bin(FILTER_CELL, & + p % coord % cell, index_tally) + if (bins(FILTER_CELL) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_CELLBORN) + ! determine next cellborn bin + bins(FILTER_CELLBORN) = get_next_bin(FILTER_CELLBORN, & + p % cell_born, index_tally) + if (bins(FILTER_CELLBORN) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_SURFACE) + ! determine next surface bin + bins(FILTER_SURFACE) = get_next_bin(FILTER_SURFACE, & + p % surface, index_tally) + if (bins(FILTER_SURFACE) == NO_BIN_FOUND) then + found_bin = .false. + return + end if + + case (FILTER_ENERGYIN) + ! determine incoming energy bin + k = t % n_filter_bins(FILTER_ENERGYIN) + ! check if energy of the particle is within energy bins + if (p % last_E < t % energy_in(1) .or. & + p % last_E > t % energy_in(k + 1)) then + found_bin = .false. + return + end if + + ! search to find incoming energy bin + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, k + 1, p % last_E) + + end select + end do FILTER_LOOP + + ! ========================================================================== + ! DETERMINE WHICH MESH CELLS TO SCORE TO + + ! Calculate number of surface crossings + n_cross = sum(abs(ijk1 - ijk0)) + 1 + + ! Copy particle's direction + uvw = p % coord0 % uvw + + ! Bounding coordinates + do j = 1, 3 + if (uvw(j) > 0) then + xyz_cross(j) = m % lower_left(j) + ijk0(j) * m % width(j) + else + xyz_cross(j) = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + end if + end do + + MESH_LOOP: do k = 1, n_cross + found_bin = .false. + + ! Calculate distance to each bounding surface. We need to treat special + ! case where the cosine of the angle is zero since this would result in a + ! divide-by-zero. + + if (k == n_cross) xyz_cross = xyz1 + + do j = 1, 3 + if (uvw(j) == 0) then + d(j) = INFINITY + else + d(j) = (xyz_cross(j) - xyz0(j))/uvw(j) + end if + end do + + ! Determine the closest bounding surface of the mesh cell by calculating + ! the minimum distance + + j = minloc(d, 1) + distance = d(j) + + ! Now use the minimum distance and diretion of the particle to determine + ! which surface was crossed + + if (all(ijk0 >= 1) .and. all(ijk0 <= m % dimension)) then + ijk_cross = ijk0 + found_bin = .true. + end if + + ! Increment indices and determine new crossing point + if (uvw(j) > 0) then + ijk0(j) = ijk0(j) + 1 + xyz_cross(j) = xyz_cross(j) + m % width(j) + else + ijk0(j) = ijk0(j) - 1 + xyz_cross(j) = xyz_cross(j) - m % width(j) + end if + + ! ======================================================================= + ! SCORE TO THIS MESH CELL + + if (found_bin) then + ! Calculate track-length estimate of flux + flux = p % wgt * distance + + ! Determine mesh bin + bins(FILTER_MESH) = mesh_indices_to_bin(m, ijk_cross) + + ! Determining scoring index + score_index = sum((bins - 1) * t % stride) + 1 + + ! Determine score for each bin + SCORE_LOOP: do i = 1, t % n_score_bins + ! determine what type of score bin + score_bin = t % score_bins(i) % scalar + + ! Determine cross section + select case(score_bin) + case (SCORE_FLUX) + score = flux + case (SCORE_TOTAL) + score = material_xs % total * flux + case (SCORE_SCATTER) + score = (material_xs % total - material_xs % absorption) * flux + case (SCORE_ABSORPTION) + score = material_xs % absorption * flux + case (SCORE_FISSION) + score = material_xs % fission * flux + case (SCORE_NU_FISSION) + score = material_xs % nu_fission * flux + end select + + ! Add score to tally + call add_to_score(t % scores(score_index, i), score) + end do SCORE_LOOP + end if + + ! Calculate new coordinates + xyz0 = xyz0 + distance * uvw + + end do MESH_LOOP + + end subroutine score_tl_on_mesh + !=============================================================================== ! GET_SCORING_BINS determines a combination of filter bins that should be scored ! for a tally based on the particle's current attributes. From 5524fde7f0d90bd7d11a29da6a7c60ebbd577102 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 25 Mar 2012 20:55:14 -0400 Subject: [PATCH 38/42] Fixed error in mesh_intersects. Closes #75 on github. --- src/mesh.F90 | 64 +++++++++++++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 3e489903b4..37847d58ac 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -218,7 +218,9 @@ contains end subroutine count_fission_sites !=============================================================================== -! MESH_INTERSECT +! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer +! boundary of the given mesh. This is important for determining whether a track +! will score to a mesh tally. !=============================================================================== function mesh_intersects(m, xyz0, xyz1) result(intersects) @@ -228,9 +230,11 @@ contains real(8), intent(in) :: xyz1(3) logical :: intersects - real(8) :: x0, y0, z0 ! track start point - real(8) :: x1, y1, z1 ! track end point - real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: x0, y0, z0 ! track start point + real(8) :: x1, y1, z1 ! track end point + real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh + real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh ! Copy coordinates of starting point x0 = xyz0(1) @@ -242,62 +246,66 @@ contains y1 = xyz1(2) z1 = xyz1(3) + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + ym0 = m % lower_left(2) + zm0 = m % lower_left(3) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + ym1 = m % upper_right(2) + zm1 = m % upper_right(3) + ! Check if line intersects bottom surface -- calculate the intersection ! point (y,z) - xi = m % lower_left(1) - yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) - if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects left surface -- calculate the intersection point ! (x,z) - yi = m % lower_left(2) - xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) - if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects front surface -- calculate the intersection point ! (x,y) - zi = m % lower_left(3) - xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) - if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then intersects = .true. return end if ! Check if line intersects top surface -- calculate the intersection ! point (y,z) - xi = m % upper_right(1) - yi = y0 + (xi - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xi - x0) * (z1 - z0) / (x1 - x0) - if (yi >= y0 .and. yi < y1 .and. zi >= z0 .and. zi < z1) then + yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects right surface -- calculate the intersection point ! (x,z) - yi = m % upper_right(2) - xi = x0 + (yi - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (yi - y0) * (z1 - z0) / (y1 - y0) - if (xi >= x0 .and. xi < x1 .and. zi >= z0 .and. zi < z1) then + xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then intersects = .true. return end if ! Check if line intersects back surface -- calculate the intersection point ! (x,y) - zi = m % upper_right(3) - xi = x0 + (zi - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zi - z0) * (y1 - y0) / (z1 - z0) - if (xi >= x0 .and. xi < x1 .and. yi >= y0 .and. yi < y1) then + xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then intersects = .true. return end if From a71840707448ec3c0521c2d709c48ce4c98b2532 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 26 Mar 2012 09:19:35 -0400 Subject: [PATCH 39/42] Fixed bug for nuclides appearing with and without S(a,b) tables. Closes #82 on github. --- src/cross_section.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b43283e037..aa155eefe7 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -188,7 +188,7 @@ contains end if ! Set last evaluated energy - micro_xs(index_nuclide) % last_E = p % E + if (index_sab == 0) micro_xs(index_nuclide) % last_E = p % E end subroutine calculate_nuclide_xs From b2c40ee3c03463bac697c88db8d072f42bdbdfd5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 26 Mar 2012 14:46:51 -0400 Subject: [PATCH 40/42] Fixed bug in incoming energy filter for track-length tallies. --- src/tally.F90 | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 1f89af152e..cb701e7070 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -727,14 +727,14 @@ contains ! determine incoming energy bin k = t % n_filter_bins(FILTER_ENERGYIN) ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(k + 1)) then + if (p % E < t % energy_in(1) .or. & + p % E > t % energy_in(k + 1)) then found_bin = .false. return end if ! search to find incoming energy bin - bins(FILTER_ENERGYIN) = binary_search(t % energy_in, k + 1, p % last_E) + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, k + 1, p % E) end select end do FILTER_LOOP @@ -857,6 +857,7 @@ contains integer :: i ! loop index for filters integer :: n ! number of bins for single filter integer :: mesh_bin ! index for mesh bin + real(8) :: E ! particle energy type(TallyObject), pointer :: t => null() type(StructuredMesh), pointer :: m => null() @@ -928,15 +929,22 @@ contains case (FILTER_ENERGYIN) ! determine incoming energy bin n = t % n_filter_bins(FILTER_ENERGYIN) + + ! make sure the correct energy is used + if (t % estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if + ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(n + 1)) then + if (E < t % energy_in(1) .or. E > t % energy_in(n + 1)) then found_bin = .false. return end if ! search to find incoming energy bin - bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + bins(FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, E) case (FILTER_ENERGYOUT) ! determine outgoing energy bin From 1466eca5fc86deae7d3d6cd9629de0cff6ad09db Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Mon, 26 Mar 2012 16:08:27 -0700 Subject: [PATCH 41/42] fixed a bug where p % last_xyz wasn't updated when necessary for tracklength tallies with meshfilters --- src/physics.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index a36202ea72..03ca45354d 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -81,6 +81,9 @@ contains ! Select smaller of the two distances distance = min(d_boundary, d_collision) + ! Save coordinates for tallying purposes + p % last_xyz = p % coord0 % xyz + ! Advance particle coord => p % coord0 do while (associated(coord)) @@ -128,9 +131,6 @@ contains p % surface = NONE call collision() - ! Save coordinates at collision for tallying purposes - p % last_xyz = p % coord0 % xyz - ! Set last material to none since cross sections will need to be ! re-evaluated p % last_material = NONE From 54feb3fe20a86edd18fc2436e4fdbaa7222925fc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 26 Mar 2012 19:51:37 -0400 Subject: [PATCH 42/42] Another fix for mesh tracklength tallies. --- src/physics.F90 | 6 +++--- src/tally.F90 | 13 +++++++++---- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 03ca45354d..ca9d389ec1 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -81,9 +81,6 @@ contains ! Select smaller of the two distances distance = min(d_boundary, d_collision) - ! Save coordinates for tallying purposes - p % last_xyz = p % coord0 % xyz - ! Advance particle coord => p % coord0 do while (associated(coord)) @@ -131,6 +128,9 @@ contains p % surface = NONE call collision() + ! Save coordinates for tallying purposes + p % last_xyz = p % coord0 % xyz + ! Set last material to none since cross sections will need to be ! re-evaluated p % last_material = NONE diff --git a/src/tally.F90 b/src/tally.F90 index cb701e7070..5d304e4ddb 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -552,7 +552,7 @@ contains ! since multiple bins can be scored to with a single track if (t % n_filter_bins(FILTER_MESH) > 0) then - call score_tl_on_mesh(tracklength_tallies(i)) + call score_tl_on_mesh(tracklength_tallies(i), distance) cycle end if @@ -621,9 +621,10 @@ contains ! these tallies, it is possible to score to multiple mesh cells for each track. !=============================================================================== - subroutine score_tl_on_mesh(index_tally) + subroutine score_tl_on_mesh(index_tally, d_track) integer, intent(in) :: index_tally + real(8), intent(in) :: d_track integer :: i ! loop index for filter/score bins integer :: j ! loop index for direction @@ -657,8 +658,8 @@ contains ! CHECK IF THIS TRACK INTERSECTS THE MESH ! Copy starting and ending location of particle - xyz0 = p % last_xyz - xyz1 = p % coord0 % xyz + xyz0 = p % coord0 % xyz - (d_track - TINY_BIT) * p % coord0 % uvw + xyz1 = p % coord0 % xyz - TINY_BIT * p % coord0 % uvw ! Determine indices for starting and ending location m => meshes(t % mesh) @@ -671,6 +672,10 @@ contains if (.not. mesh_intersects(m, xyz0, xyz1)) return end if + ! Reset starting and ending location + xyz0 = p % coord0 % xyz - d_track * p % coord0 % uvw + xyz1 = p % coord0 % xyz + ! ========================================================================= ! CHECK FOR SCORING COMBINATION FOR FILTERS OTHER THAN MESH