diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst index 874290465d..7105e9f7d1 100644 --- a/docs/source/usersguide/processing.rst +++ b/docs/source/usersguide/processing.rst @@ -288,7 +288,7 @@ SILO file can be created with: import silomesh as sm sm.init_silo("fluxtally.silo") - sm.init_mesh('tally_mesh',*mesh.dimension, *mesh.lower_left, *mesh.upper_right) + sm.init_mesh('tally_mesh', *mesh.dimension, *mesh.lower_left, *mesh.upper_right) sm.init_var('flux_tally_thermal') for x in range(1,nx+1): for y in range(1,ny+1): diff --git a/docs/source/usersguide/troubleshoot.rst b/docs/source/usersguide/troubleshoot.rst index 64b955593e..49cf41c8de 100644 --- a/docs/source/usersguide/troubleshoot.rst +++ b/docs/source/usersguide/troubleshoot.rst @@ -98,8 +98,38 @@ with the :envvar:`CROSS_SECTIONS` environment variable. It is recommended to add a line in your ``.profile`` or ``.bash_profile`` setting the :envvar:`CROSS_SECTIONS` environment variable. +Geometry Debugging +****************** + +Overlapping Cells +^^^^^^^^^^^^^^^^^ + +For fast run times, normal simulations do not check if the geometry is +incorrectly defined to have overlapping cells. This can lead to incorrect +results that may or may not be obvious when there are errors in the geometry +input file. The built-in 2D and 3D plotters will check for cell overlaps at +the center of every pixel or voxel position they process, however this might +not be a sufficient check to ensure correctly defined geometry. For instance, +if an overlap is of small aspect ratio, the plotting resolution might not be +high enough to produce any pixels in the overlapping area. + +To reliably validate a geometry input, it is best to run the problem in +geometry debugging mode with the ``-g``, ``-geometry-debug``, or +``--geometry-debug`` command-line options. This will enable checks for +overlapping cells at every move of esch simulated particle. Depending on the +complexity of the geometry input file, this could add considerable overhead to +the run (these runs can still be done in parallel). As a result, for this run +mode the user will probably want to run fewer particles than a normal +simulation run. In this case it is important to be aware of how much coverage +each area of the geometry is getting. For instance, if certain regions do not +have many particles travelling through them there will not be many locations +where overlaps are checked for in that region. The user should refer to the +output after a geometry debug run to see how many checks were performed in each +cell, and then adjust the number of starting particles or starting source +distributions accordingly to achieve good coverage. + ERROR: After particle __ crossed surface __ it could not be located in any cell and it did not leak. -**************************************************************************************************** +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This error can arise either if a problem is specified with no boundary conditions or if there is an error in the geometry itself. First check to ensure @@ -120,5 +150,10 @@ has a collision. For example, if you received this error at cycle 5, generation 5 1 4032 +For large runs it is often advantageous to run only the offending particle by +using particle restart mode with the ``-s``, ``-particle``, or ``--particle`` +command-line options in conjunction with the particle restart files that are +created when particles are lost with this error. + .. _gfortran: http://gcc.gnu.org/wiki/GFortran .. _mailing list: https://groups.google.com/forum/?fromgroups=#!forum/openmc-users diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index b8a9d9c44c..aefbca9171 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -13,6 +13,10 @@ It is assumed that if no is specified, the XML input files are present in the current directory. .SH OPTIONS .TP +.B "\-g\fR, \fP\-\-geometry-debug" +Run in geometry debugging mode, where cell overlaps are checked for after each +move of a particle +.TP .B "\-p\fR, \fP\-\-plot" Run in plotting mode .TP diff --git a/src/finalize.F90 b/src/finalize.F90 index bf0c752ffb..adffe7b711 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -4,7 +4,8 @@ module finalize use cmfd_output, only: finalize_cmfd # endif use global - use output, only: print_runtime, print_results, write_tallies + use output, only: print_runtime, print_results, & + print_overlap_check, write_tallies use tally, only: tally_statistics #ifdef MPI @@ -35,6 +36,7 @@ contains if (output_tallies) then if (master) call write_tallies() end if + if (check_overlaps) call reduce_overlap_count() end if #ifdef PETSC @@ -50,6 +52,7 @@ contains run_mode /= MODE_PARTICLE)) then call print_runtime() call print_results() + if (check_overlaps) call print_overlap_check() end if ! deallocate arrays @@ -71,4 +74,22 @@ contains end subroutine finalize_run +!=============================================================================== +! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master +!=============================================================================== + + subroutine reduce_overlap_count() + +#ifdef MPI + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, & + MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, & + MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + end if +#endif + + end subroutine reduce_overlap_count + end module finalize diff --git a/src/geometry.F90 b/src/geometry.F90 index 448691ca2d..cb90997846 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -62,6 +62,58 @@ contains end function simple_cell_contains + +!=============================================================================== +! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's +! position using simple_cell_contains and the LocalCoord's built up by find_cell +!=============================================================================== + + subroutine check_cell_overlap() + + integer :: i ! cell loop index on a level + integer :: n ! number of cells to search on a level + integer :: index_cell ! index in cells array + type(Cell), pointer :: c ! pointer to cell + type(Universe), pointer :: univ ! universe to search in + type(LocalCoord), pointer :: coord ! particle coordinate to search on + + coord => p % coord0 + + ! loop through each coordinate level + do while (associated(coord)) + + p % coord => coord + + univ => universes(coord % universe) + n = univ % n_cells + + ! loop through each cell on this level + do i = 1, n + index_cell = univ % cells(i) + c => cells(index_cell) + + if (simple_cell_contains(c)) then + ! the particle should only be contained in one cell per level + if (index_cell /= coord % cell) then + message = "Overlapping cells detected: " // & + trim(to_str(cells(index_cell) % id)) // ", " // & + trim(to_str(cells(coord % cell) % id)) // & + " on universe " // trim(to_str(univ % id)) + call fatal_error() + end if + + overlap_check_cnt(index_cell) = overlap_check_cnt(index_cell) + 1 + + end if + + end do + + coord => coord % next + + end do + + end subroutine check_cell_overlap + !=============================================================================== ! FIND_CELL determines what cell a source particle is in within a particular ! universe. If the base universe is passed, the particle should be found as long diff --git a/src/global.F90 b/src/global.F90 index 5597ae2dfb..ee9f752295 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -274,6 +274,10 @@ module global ! screen and in logs integer :: verbosity = 7 + ! Flag for enabling cell overlap checking during transport + logical :: check_overlaps = .false. + integer(8), allocatable :: overlap_check_cnt(:) + ! Trace for single particle logical :: trace integer :: trace_batch @@ -396,6 +400,9 @@ contains if (allocated(materials)) deallocate(materials) if (allocated(plots)) deallocate(plots) + ! Deallocate geometry debugging information + if (allocated(overlap_check_cnt)) deallocate(overlap_check_cnt) + ! Deallocate cross section data, listings, and cache if (allocated(nuclides)) then ! First call the clear routines diff --git a/src/initialize.F90 b/src/initialize.F90 index a04c8a2200..2216a4c4b1 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -5,14 +5,15 @@ module initialize use constants use dict_header, only: DictIntInt, ElemKeyValueII use energy_grid, only: unionized_grid - use error, only: fatal_error + use error, only: fatal_error, warning use geometry, only: neighbor_lists 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 use output, only: title, header, write_summary, print_version, & - print_usage, write_xs_summary, print_plot + print_usage, write_xs_summary, print_plot, & + write_message use output_interface, only: file_open, file_close, read_data use random_lcg, only: initialize_prng use source, only: initialize_source @@ -142,9 +143,17 @@ contains end if end if - ! check for particle restart run + ! Check for particle restart run if (particle_restart_run) run_mode = MODE_PARTICLE + ! Warn if overlap checking is on + if (master .and. check_overlaps) then + message = "" + call write_message() + message = "Cell overlap checking is ON" + call warning() + end if + ! Stop initialization timer call time_initialize % stop() @@ -316,6 +325,8 @@ contains select case (argv(i)) case ('-p', '-plot', '--plot') run_mode = MODE_PLOTTING + check_overlaps = .true. + case ('-n', '-n_particles', '--n_particles') ! Read number of particles per cycle i = i + 1 @@ -354,6 +365,9 @@ contains path_state_point = argv(i) restart_run = .true. + case ('-g', '-geometry-debug', '--geometry-debug') + check_overlaps = .true. + case ('-?', '-help', '--help') call print_usage() stop diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a9b62d4905..45abaf9a2e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -626,6 +626,11 @@ contains ! Allocate cells array allocate(cells(n_cells)) + if (check_overlaps) then + allocate(overlap_check_cnt(n_cells)) + overlap_check_cnt = 0 + end if + n_universes = 0 do i = 1, n_cells c => cells(i) diff --git a/src/output.F90 b/src/output.F90 index b8302226d7..b3481f509b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -165,6 +165,7 @@ contains write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]' write(OUTPUT_UNIT,*) write(OUTPUT_UNIT,*) 'Options:' + write(OUTPUT_UNIT,*) ' -g, --geometry-debug Run in geometry debugging mode' write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point' write(OUTPUT_UNIT,*) ' or a particle restart file' @@ -1478,6 +1479,42 @@ contains end subroutine print_results +!=============================================================================== +! PRINT_OVERLAP_DEBUG displays information regarding overlap checking results +!=============================================================================== + + subroutine print_overlap_check + + integer :: i, j + integer :: num_sparse = 0 + + ! display header block for geometry debugging section + call header("Cell Overlap Check Summary") + + write(ou,100) 'Cell ID','No. Overlap Checks' + + do i = 1, n_cells + write(ou,101) cells(i) % id, overlap_check_cnt(i) + if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1 + end do + write(ou,*) + write(ou,'(1X,A)') 'There were ' // trim(to_str(num_sparse)) // & + ' cells with less than 10 overlap checks' + j = 0 + do i = 1, n_cells + if (overlap_check_cnt(i) < 10) then + j = j + 1 + write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id)) + if (modulo(j,8) == 0) write(ou,*) + end if + end do + write(ou,*) + +100 format (1X,A,T15,A) +101 format (1X,I8,T15,I12) + + end subroutine print_overlap_check + !=============================================================================== ! WRITE_TALLIES creates an output file and writes out the mean values of all ! tallies and their standard deviations diff --git a/src/physics.F90 b/src/physics.F90 index 707049b1a8..4dc5e13b29 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -7,7 +7,8 @@ module physics use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed use geometry, only: find_cell, distance_to_boundary, & - cross_surface, cross_lattice + cross_surface, cross_lattice, & + check_cell_overlap use geometry_header, only: Universe, BASE_UNIVERSE use global use interpolation, only: interpolate_tab1 @@ -75,6 +76,8 @@ contains do while (p % alive) + if (check_overlaps) call check_cell_overlap() + ! Calculate microscopic and macroscopic cross sections -- note: if the ! material is the same as the last material and the energy of the ! particle hasn't changed, we don't need to lookup cross sections again. diff --git a/src/plot.F90 b/src/plot.F90 index c38788402a..e7b8a235a6 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -2,8 +2,8 @@ module plot use constants use error, only: fatal_error - use geometry, only: find_cell - use geometry_header, only: Universe, BASE_UNIVERSE + use geometry, only: find_cell, check_cell_overlap + use geometry_header, only: Cell, BASE_UNIVERSE use global use output, only: write_message use particle_header, only: deallocate_coord @@ -62,6 +62,7 @@ contains p % coord => p % coord0 call find_cell(found_cell) + if (check_overlaps) call check_cell_overlap() if (.not. found_cell) then ! If no cell, revert to default color