Merge pull request #573 from paulromano/subcritical

Allow fission in fixed source simulations

src/particle_header.F90

@@ -86,7 +86,7 @@ module particle_header
     logical    :: write_track = .false.
 
     ! Secondary particles created
-    integer    :: n_secondary = 0
+    integer(8) :: n_secondary = 0
     type(Bank) :: secondary_bank(MAX_SECONDARY)
 
   contains

src/physics.F90

@@ -88,9 +88,14 @@ subroutine sample_reaction(p)
     ! change when sampling fission sites. The following block handles all
     ! absorption (including fission)
 
-    if (nuc % fissionable .and. run_mode == MODE_EIGENVALUE) then
+    if (nuc % fissionable) then
       call sample_fission(i_nuclide, i_reaction)
-      call create_fission_sites(p, i_nuclide, i_reaction)
+      if (run_mode == MODE_EIGENVALUE) then
+        call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank)
+      elseif (run_mode == MODE_FIXEDSOURCE) then
+        call create_fission_sites(p, i_nuclide, i_reaction, &
+             p % secondary_bank, p % n_secondary)
+      end if
     end if
 
     ! If survival biasing is being used, the following subroutine adjusts the
@@ -1071,10 +1076,12 @@ end subroutine sample_cxs_target_velocity
 ! neutrons produced from fission and creates appropriate bank sites.
 !===============================================================================
 
-  subroutine create_fission_sites(p, i_nuclide, i_reaction)
+  subroutine create_fission_sites(p, i_nuclide, i_reaction, bank_array, size_bank)
     type(Particle), intent(inout) :: p
     integer,        intent(in)    :: i_nuclide
     integer,        intent(in)    :: i_reaction
+    type(Bank),     intent(inout) :: bank_array(:)
+    integer(8),     intent(inout) :: size_bank
 
     integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
     integer :: i                        ! loop index
@@ -1123,27 +1130,35 @@ subroutine create_fission_sites(p, i_nuclide, i_reaction)
       nu = int(nu_t) + 1
     end if
 
-    ! Check for fission bank size getting hit
-    if (n_bank + nu > size(fission_bank)) then
-      if (master) call warning("Maximum number of sites in fission bank &
-           &reached. This can result in irreproducible results using different &
-           &numbers of processes/threads.")
+    ! Check for bank size getting hit. For fixed source calculations, this is a
+    ! fatal error. For eigenvalue calculations, it just means that k-effective
+    ! was too high for a single batch.
+    if (size_bank + nu > size(bank_array)) then
+      if (run_mode == MODE_FIXEDSOURCE) then
+        call fatal_error("Secondary particle bank size limit reached. If you &
+             &are running a subcritical multiplication problem, k-effective &
+             &may be too close to one.")
+      else
+        if (master) call warning("Maximum number of sites in fission bank &
+             &reached. This can result in irreproducible results using different &
+             &numbers of processes/threads.")
+      end if
     end if
 
     ! Bank source neutrons
-    if (nu == 0 .or. n_bank == size(fission_bank)) return
+    if (nu == 0 .or. size_bank == size(bank_array)) return
 
     ! Initialize counter of delayed neutrons encountered for each delayed group
     ! to zero.
     nu_d(:) = 0
 
     p % fission = .true. ! Fission neutrons will be banked
-    do i = int(n_bank,4) + 1, int(min(n_bank + nu, int(size(fission_bank),8)),4)
+    do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
       ! Bank source neutrons by copying particle data
-      fission_bank(i) % xyz = p % coord(1) % xyz
+      bank_array(i) % xyz = p % coord(1) % xyz
 
       ! Set weight of fission bank site
-      fission_bank(i) % wgt = ONE/weight
+      bank_array(i) % wgt = ONE/weight
 
       ! Sample cosine of angle -- fission neutrons are always emitted
       ! isotropically. Sometimes in ACE data, fission reactions actually have
@@ -1153,17 +1168,17 @@ subroutine create_fission_sites(p, i_nuclide, i_reaction)
 
       ! Sample azimuthal angle uniformly in [0,2*pi)
       phi = TWO*PI*prn()
-      fission_bank(i) % uvw(1) = mu
-      fission_bank(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
-      fission_bank(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
+      bank_array(i) % uvw(1) = mu
+      bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
+      bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
 
       ! Sample secondary energy distribution for fission reaction and set energy
       ! in fission bank
-      fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(&
-           i_reaction), p)
+      bank_array(i) % E = sample_fission_energy(nuc, &
+           nuc % reactions(i_reaction), p)
 
       ! Set the delayed group of the neutron
-      fission_bank(i) % delayed_group = p % delayed_group
+      bank_array(i) % delayed_group = p % delayed_group
 
       ! Increment the number of neutrons born delayed
       if (p % delayed_group > 0) then
@@ -1172,7 +1187,7 @@ subroutine create_fission_sites(p, i_nuclide, i_reaction)
     end do
 
     ! increment number of bank sites
-    n_bank = min(n_bank + nu, int(size(fission_bank),8))
+    size_bank = min(size_bank + nu, int(size(bank_array),8))
 
     ! Store total and delayed weight banked for analog fission tallies
     p % n_bank   = nu

src/tracking.F90

@@ -202,6 +202,7 @@ subroutine transport(p)
         if (p % n_secondary > 0) then
           call p % initialize_from_source(p % secondary_bank(p % n_secondary))
           p % n_secondary = p % n_secondary - 1
+          n_event = 0
 
           ! Enter new particle in particle track file
           if (p % write_track) call add_particle_track()

tests/test_fixed_source/materials.xml

@@ -4,6 +4,7 @@
   <material id="1">
     <density value="7.5" units="g/cc" />
     <nuclide name="O-16" xs="71c" ao="1.0" />
+    <nuclide name="U-238" xs="71c" ao="0.0001" />
   </material>
 
 </materials>

tests/test_fixed_source/results_true.dat

@@ -1,6 +1,6 @@
 tally 1:
-4.538791E+02
-2.073271E+04
+4.563929E+02
+2.091711E+04
 leakage:
-9.830000E+00
-9.663900E+00
+9.780000E+00
+9.566400E+00