Merge pull request #312 from bhermanmit/fission_src

Added capability to filter box source by fissionable material

docs/source/usersguide/input.rst

@@ -311,6 +311,12 @@ attributes/sub-elements:
       parallelepiped and the last three of which specify the upper-right
       corner. Source sites are sampled uniformly through that parallelepiped.
 
+      To filter a "box" spatial distribution by fissionable material, specify
+      "fission" tag instead of "box". The ``parameters`` should be given as six
+      real numbers, the first three of which specify the lower-left corner of a
+      parallelepiped and the last three of which specify the upper-right
+      corner. Source sites are sampled uniformly through that parallelepiped.
+
       For a "point" spatial distribution, ``parameters`` should be given as
       three real numbers which specify the (x,y,z) location of an isotropic
       point source

src/ace.F90

@@ -172,6 +172,25 @@ contains
     ! Avoid some valgrind leak errors
     call already_read % clear()
 
+    ! Loop around material
+    MATERIAL_LOOP3: do i = 1, n_materials
+
+      ! Get material
+      mat => materials(i)
+
+      ! Loop around nuclides in material
+      NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
+
+        ! Check for fission in nuclide
+        if (nuclides(mat % nuclide(j)) % fissionable) then
+          mat % fissionable = .true.
+          exit NUCLIDE_LOOP2
+        end if
+
+      end do NUCLIDE_LOOP2
+
+    end do MATERIAL_LOOP3
+
   end subroutine read_xs
 
 !===============================================================================

src/constants.F90

@@ -335,8 +335,9 @@ module constants
 
   ! Source spatial distribution types
   integer, parameter :: &
-       SRC_SPACE_BOX   = 1, & ! Source in a rectangular prism
-       SRC_SPACE_POINT = 2    ! Source at a single point
+       SRC_SPACE_BOX     = 1, & ! Source in a rectangular prism
+       SRC_SPACE_POINT   = 2, & ! Source at a single point
+       SRC_SPACE_FISSION = 3    ! Source in prism filtered by fissionable mats
 
   ! Source angular distribution types
   integer, parameter :: &

src/initialize.F90

@@ -11,6 +11,7 @@ module initialize
   use global
   use input_xml,        only: read_input_xml, read_cross_sections_xml,         &
                               cells_in_univ_dict, read_plots_xml
+  use material_header,  only: Material
   use output,           only: title, header, write_summary, print_version,     &
                               print_usage, write_xs_summary, print_plot,       &
                               write_message

src/input_xml.F90

@@ -292,6 +292,9 @@ contains
         case ('box')
           external_source % type_space = SRC_SPACE_BOX
           coeffs_reqd = 6
+        case ('fission')
+          external_source % type_space = SRC_SPACE_FISSION
+          coeffs_reqd = 6
         case ('point')
           external_source % type_space = SRC_SPACE_POINT
           coeffs_reqd = 3

src/material_header.F90

@@ -21,6 +21,10 @@ module material_header
     ! Temporary names read during initialization
     character(12), allocatable :: names(:)     ! isotope names
     character(12), allocatable :: sab_names(:) ! name of S(a,b) table
+
+    ! Does this material contain fissionable nuclides?
+    logical :: fissionable = .false.
+
   end type Material
 
 end module material_header

src/source.F90

@@ -135,6 +135,42 @@ contains
       end do
       call p % clear()
 
+    case (SRC_SPACE_FISSION)
+      ! Repeat sampling source location until a good site has been found
+      found = .false.
+      do while (.not.found)
+        ! Set particle defaults
+        call p % initialize()
+
+        ! 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
+        p % coord0 % uvw = [ ONE, ZERO, ZERO ]
+
+        ! 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
+          cycle
+        end if
+        if (p % material == MATERIAL_VOID) then
+          found = .false.
+          cycle
+        end if
+        if (.not. materials(p % material) % fissionable) found = .false.
+      end do
+      call p % clear()
+
     case (SRC_SPACE_POINT)
       ! Point source
       site % xyz = external_source % params_space