Merge pull request #110 from valeriaRaffuzzi/varRed

Adding MGimp processor and maximum number of split per particle
CambridgeNuclear · Jan 25, 2024 · 16ee2a7 · 16ee2a7
2 parents 58df62d + 82de4c7
commit 16ee2a7
Show file tree

Hide file tree

Showing 8 changed files with 450 additions and 15 deletions.
diff --git a/CollisionOperator/CollisionProcessors/CMakeLists.txt b/CollisionOperator/CollisionProcessors/CMakeLists.txt
@@ -2,5 +2,6 @@
 add_sources( ./collisionProcessor_inter.f90 
              ./collisionProcessorFactory_func.f90 
              ./neutronCEstd_class.f90
-	     ./neutronCEimp_class.f90
-             ./neutronMGstd_class.f90) 
+             ./neutronCEimp_class.f90
+             ./neutronMGstd_class.f90
+             ./neutronMGimp_class.f90)
diff --git a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90
@@ -11,6 +11,7 @@ module collisionProcessorFactory_func
   use neutronCEstd_class, only : neutronCEstd
   use neutronCEimp_class, only : neutronCEimp
   use neutronMGstd_class, only : neutronMGstd
+  use neutronMGimp_class, only : neutronMGimp
 
   implicit none
   private
@@ -23,7 +24,8 @@ module collisionProcessorFactory_func
   ! For now  it is necessary to adjust trailing blanks so all enteries have the same length
   character(nameLen),dimension(*),parameter :: AVALIBLE_collisionProcessors = [ 'neutronCEstd',&
                                                                                 'neutronCEimp',&
-                                                                                'neutronMGstd']
+                                                                                'neutronMGstd',&
+                                                                                'neutronMGimp']
 
 contains
 
@@ -54,6 +56,10 @@ subroutine new_collisionProcessor(new,dict)
       case('neutronMGstd')
         allocate(neutronMGstd :: new)
 
+      case('neutronMGimp')
+        allocate(neutronMGimp :: new)
+        call new % init(dict)
+
       case default
         print *, AVALIBLE_collisionProcessors
         call fatalError(Here, 'Unrecognised type of collisionProcessor: ' // trim(type))

diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90
@@ -122,7 +122,9 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
     ! Note: the ordering must not be changed between feeding the particle to the tally
     ! and updating the particle's preCollision state, otherwise this may cause certain
     ! tallies (e.g., collisionProbability) to return dubious results
+
     call tally % reportInColl(p, .false.)
+
     call p % savePreCollision()
 
     ! Choose collision nuclide and general type (Scatter, Capture or Fission)

diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
@@ -62,6 +62,7 @@ module neutronCEimp_class
   !!  impAbs  -> is implicit capture performed? (off by default)
   !!  impGen  -> are fission sites generated implicitly? (on by default)
   !!  UFS     -> uniform fission sites variance reduction
+  !!  maxSplit -> maximum number of splits allowed per particle (default = 1000)
   !!  thresh_E -> Energy threshold for explicit treatment of target nuclide movement [-].
   !!              Target movment is sampled if neutron energy E < kT * thresh_E where
   !!              kT is target material temperature in [MeV]. (default = 400.0)
@@ -89,6 +90,7 @@ module neutronCEimp_class
   !!   #impGen         <logical>;#
   !!   #UFS            <logical>;#
   !!   #weightWindows  <logical>;#
+  !!   #maxSplit       <integer>;#
   !!   }
   !!
   type, public, extends(collisionProcessor) :: neutronCEimp
@@ -109,14 +111,15 @@ module neutronCEimp_class
     real(defReal) :: thresh_A
     real(defReal) :: DBRCeMin
     real(defReal) :: DBRCeMax
+    integer(shortInt) :: maxSplit
 
     ! Variance reduction options
-    logical(defBool) :: weightWindows
-    logical(defBool) :: splitting
-    logical(defBool) :: roulette
-    logical(defBool) :: implicitAbsorption ! Prevents particles dying through capture
-    logical(defBool) :: implicitSites ! Generates fission sites on every fissile collision
-    logical(defBool) :: uniFissSites
+    logical(defBool)  :: weightWindows
+    logical(defBool)  :: splitting
+    logical(defBool)  :: roulette
+    logical(defBool)  :: implicitAbsorption ! Prevents particles dying through capture
+    logical(defBool)  :: implicitSites ! Generates fission sites on every fissile collision
+    logical(defBool)  :: uniFissSites
 
     type(weightWindowsField), pointer :: weightWindowsMap
 
@@ -168,6 +171,7 @@ subroutine init(self, dict)
 
     ! Obtain settings for variance reduction
     call dict % getOrDefault(self % weightWindows,'weightWindows', .false.)
+    call dict % getOrDefault(self % maxSplit,'maxSplit', 1000)
     call dict % getOrDefault(self % splitting,'split', .false.)
     call dict % getOrDefault(self % roulette,'roulette', .false.)
     call dict % getOrDefault(self % minWgt,'minWgt',0.25_defReal)
@@ -194,8 +198,8 @@ subroutine init(self, dict)
     end if
 
     if (self % implicitAbsorption) then
-      if (.not.self % roulette) call fatalError(Here,&
-         'Must use Russian roulette when using implicit absorption')
+      if (.not.self % roulette .and. .not. self % weightWindows) call fatalError(Here,&
+         'Must use Russian roulette or weight windows when using implicit absorption')
       if (.not.self % implicitSites) call fatalError(Here,&
          'Must generate fission sites implicitly when using implicit absorption')
     end if
@@ -533,7 +537,8 @@ subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle)
     real(defReal)                        :: minWgt, maxWgt, avWgt
 
     if (p % isDead) then
-    ! Do nothing
+
+      ! Do nothing !
 
     elseif (p % E < self % minE) then
       p % isDead = .true.
@@ -547,10 +552,13 @@ subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle)
 
       ! If a particle is outside the WW map and all the weight limits
       ! are zero nothing happens. NOTE: this holds for positive weights only
-      if ((p % w > maxWgt) .and. (maxWgt /= ZERO)) then
+
+      if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < self % maxSplit)) then
         call self % split(p, tally, thisCycle, maxWgt)
+
       elseif (p % w < minWgt) then
         call self % russianRoulette(p, avWgt)
+
       end if
 
     ! Splitting with fixed threshold
@@ -596,7 +604,14 @@ subroutine split(self, p, tally, thisCycle, maxWgt)
     ! This value must be at least 2
     mult = ceiling(p % w/maxWgt)
 
+    ! Limit maximum split
+    if (mult + p % splitCount > self % maxSplit) then
+      mult = self % maxSplit - p % splitCount + 1
+    end if
+
     ! Copy particle to a particle state
+    ! Note that particleState doesn't have property splitCount, so it is reset
+    ! to 0 for the new particle
     pTemp = p
     pTemp % wgt = p % w/mult
 
@@ -606,6 +621,9 @@ subroutine split(self, p, tally, thisCycle, maxWgt)
       call tally % reportSpawn(N_N_SPLIT, p, pTemp)
     end do
 
+    ! Update particle split count
+    p % splitCount = p % splitCount + mult
+
     ! Decrease original particle weight
     p % w = p % w/mult