Merge branch 'cell_overlap_check'

Conflicts:
	src/initialize.F90

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):

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
 
     <trace>5 1 4032</trace>
 
+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

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

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

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

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

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

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)

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

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.

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