Merge pull request #225 from nelsonag/resample

Resampling for starting source (to close issue #73)

src/DEPENDENCIES

@@ -329,6 +329,7 @@ solver_interface.o: vector_header.o
 source.o: bank_header.o
 source.o: constants.o
 source.o: error.o
+source.o: geometry.o
 source.o: geometry_header.o
 source.o: global.o
 source.o: math.o

src/constants.F90

@@ -20,7 +20,7 @@ module constants
        FILETYPE_PARTICLE_RESTART = -2
 
   ! ============================================================================
-  ! ADJUSTABLE PARAMETERS 
+  ! ADJUSTABLE PARAMETERS
 
   ! NOTE: This is the only section of the constants module that should ever be
   ! adjusted. Modifying constants in other sections may cause the code to fail.
@@ -50,6 +50,10 @@ module constants
   integer, parameter :: MAX_WORD_LEN = 150
   integer, parameter :: MAX_FILE_LEN = 255
 
+  ! Maximum number of external source spatial resamples to encounter before an
+  ! error is thrown.
+  integer, parameter :: MAX_EXTSRC_RESAMPLES = 10000
+
   ! ============================================================================
   ! PHYSICAL CONSTANTS
 
@@ -110,9 +114,9 @@ module constants
 
   ! Surface types
   integer, parameter ::  &
-       SURF_PX     =  1, & ! Plane parallel to x-plane 
-       SURF_PY     =  2, & ! Plane parallel to y-plane 
-       SURF_PZ     =  3, & ! Plane parallel to z-plane 
+       SURF_PX     =  1, & ! Plane parallel to x-plane
+       SURF_PY     =  2, & ! Plane parallel to y-plane
+       SURF_PZ     =  3, & ! Plane parallel to z-plane
        SURF_PLANE  =  4, & ! Arbitrary plane
        SURF_CYL_X  =  5, & ! Cylinder along x-axis
        SURF_CYL_Y  =  6, & ! Cylinder along y-axis
@@ -143,7 +147,7 @@ module constants
        ELECTRON = 3
 
   ! Angular distribution type
-  integer, parameter :: & 
+  integer, parameter :: &
        ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution
        ANGLE_32_EQUI   = 2, & ! 32 equiprobable bins
        ANGLE_TABULAR   = 3    ! Tabular angular distribution
@@ -275,7 +279,7 @@ module constants
        SCORE_KAPPA_FISSION = -12, & ! fission energy production rate
        SCORE_CURRENT       = -13, & ! partial current
        SCORE_EVENTS        = -14    ! number of events
-       
+
   ! Maximum scattering order supported
   integer, parameter :: SCATT_ORDER_MAX = 10
   character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10"
@@ -322,7 +326,7 @@ module constants
 
   ! Source angular distribution types
   integer, parameter :: &
-       SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular 
+       SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular
        SRC_ANGLE_MONO      = 2, & ! Monodirectional source
        SRC_ANGLE_TABULAR   = 3    ! Tabular distribution
 
@@ -332,7 +336,7 @@ module constants
        SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum
        SRC_ENERGY_WATT    = 3, & ! Watt fission spectrum
        SRC_ENERGY_TABULAR = 4    ! Tabular distribution
-       
+
   ! ============================================================================
   ! MISCELLANEOUS CONSTANTS
 

src/source.F90

@@ -3,6 +3,7 @@ module source
   use bank_header,     only: Bank
   use constants
   use error,           only: fatal_error
+  use geometry,        only: find_cell
   use geometry_header, only: BASE_UNIVERSE
   use global
   use math,            only: maxwell_spectrum, watt_spectrum
@@ -73,6 +74,9 @@ contains
     real(8) :: p_max(3)   ! maximum coordinates of source
     real(8) :: a          ! Arbitrary parameter 'a'
     real(8) :: b          ! Arbitrary parameter 'b'
+    logical :: found      ! Does the source particle exist within geometry?
+    type(Particle) :: p   ! Temporary particle for using find_cell
+    integer, save :: num_resamples = 0 ! Number of resamples encountered
 
     ! Set weight to one by default
     site % wgt = ONE
@@ -80,11 +84,31 @@ contains
     ! Sample position
     select case (external_source % type_space)
     case (SRC_SPACE_BOX)
-      ! Coordinates sampled uniformly over a box
-      p_min = external_source % params_space(1:3)
-      p_max = external_source % params_space(4:6)
-      r = (/ (prn(), i = 1,3) /)
-      site % xyz = p_min + r*(p_max - p_min)
+      ! Set particle defaults
+      call p % initialize()
+      ! Repeat sampling source location until a good site has been found
+      found = .false.
+      do while (.not.found)
+        ! Coordinates sampled uniformly over a box
+        p_min = external_source % params_space(1:3)
+        p_max = external_source % params_space(4:6)
+        r = (/ (prn(), i = 1,3) /)
+        site % xyz = p_min + r*(p_max - p_min)
+
+        ! Fill p with needed data
+        p % coord0 % xyz = site % xyz
+
+        ! Now search to see if location exists in geometry
+        call find_cell(p, found)
+        if (.not. found) then
+          num_resamples = num_resamples + 1
+          if (num_resamples == MAX_EXTSRC_RESAMPLES) then
+            message = "Maximum number of external source spatial resamples &
+                      &reached!"
+            call fatal_error()
+          end if
+        end if
+      end do
 
     case (SRC_SPACE_POINT)
       ! Point source
@@ -146,7 +170,7 @@ contains
   end subroutine sample_external_source
 
 !===============================================================================
-! GET_SOURCE_PARTICLE returns the next source particle 
+! GET_SOURCE_PARTICLE returns the next source particle
 !===============================================================================
 
   subroutine get_source_particle(p, index_source)