Merge branch 'DTmajorant' of https://github.com/valeriaRaffuzzi/SCONE …

…into DTmajorant

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)

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

CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90

@@ -86,13 +86,15 @@ module collisionProcessor_inter
     !! Procedure interface for all customisable actions associated with
     !! processing of sollision event (scatter, fission etc.)
     !!
-    subroutine collisionAction(self, p, collDat, thisCycle, nextCycle)
+    subroutine collisionAction(self, p, tally, collDat, thisCycle, nextCycle)
       import :: collisionProcessor, &
                 collisionData, &
+                tallyAdmin, &
                 particle,&
                 particleDungeon
       class(collisionProcessor), intent(inout) :: self
       class(particle), intent(inout)           :: p
+      type(tallyAdmin), intent(inout)          :: tally
       type(collisionData), intent(inout)       :: collDat
       class(particleDungeon),intent(inout)     :: thisCycle
       class(particleDungeon),intent(inout)     :: nextCycle
@@ -104,7 +106,7 @@ end subroutine collisionAction
   !!
   !! Generic flow of collision processing
   !!
-  subroutine collide(self, p, tally ,thisCycle, nextCycle)
+  subroutine collide(self, p, tally, thisCycle, nextCycle)
     class(collisionProcessor), intent(inout) :: self
     class(particle), intent(inout)           :: p
     type(tallyAdmin), intent(inout)          :: tally
@@ -118,30 +120,32 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle)
 
     ! Report in-collision & save pre-collison state
     ! 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 
+    ! 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)
-    call self % sampleCollision(p, collDat, thisCycle, nextCycle)
+    call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)
 
     ! Perform implicit treatment
-    call self % implicit(p, collDat, thisCycle, nextCycle)
+    call self % implicit(p, tally, collDat, thisCycle, nextCycle)
 
     ! Select physics to be processed based on MT number
     select case(collDat % MT)
       case(N_N_elastic, macroAllScatter)
-        call self % elastic(p, collDat, thisCycle, nextCycle)
+        call self % elastic(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_N_inelastic, macroIEScatter)
-        call self % inelastic(p, collDat, thisCycle, nextCycle)
+        call self % inelastic(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_DISAP, macroCapture)
-        call self % capture(p, collDat, thisCycle, nextCycle)
+        call self % capture(p, tally, collDat, thisCycle, nextCycle)
 
       case(N_FISSION, macroFission)
-        call self % fission(p, collDat, thisCycle, nextCycle)
+        call self % fission(p, tally, collDat, thisCycle, nextCycle)
 
       case(noInteraction)
         ! Do nothing
@@ -152,7 +156,7 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle)
     end select
 
     ! Apply post collision implicit treatments
-    call self % cutoffs(p, collDat, thisCycle, nextCycle)
+    call self % cutoffs(p, tally, collDat, thisCycle, nextCycle)
 
     ! Update particle collision counter
     p % collisionN = p % collisionN + 1

CollisionOperator/CollisionProcessors/neutronCEimp_class.f90

@@ -38,6 +38,10 @@ module neutronCEimp_class
   ! Scattering procedures
   use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, &
                                      asymptoticInelasticScatter, targetVelocity_DBRCXS
+
+  ! Tally interfaces
+  use tallyAdmin_class,       only : tallyAdmin
+
   implicit none
   private
 
@@ -58,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)
@@ -85,6 +90,7 @@ module neutronCEimp_class
   !!   #impGen         <logical>;#
   !!   #UFS            <logical>;#
   !!   #weightWindows  <logical>;#
+  !!   #maxSplit       <integer>;#
   !!   }
   !!
   type, public, extends(collisionProcessor) :: neutronCEimp
@@ -105,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
 
@@ -164,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)
@@ -190,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
@@ -213,9 +221,10 @@ end subroutine init
   !!
   !! Samples collision without any implicit treatment
   !!
-  subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
+  subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)   :: self
     class(particle), intent(inout)       :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)   :: collDat
     class(particleDungeon),intent(inout) :: thisCycle
     class(particleDungeon),intent(inout) :: nextCycle
@@ -252,9 +261,10 @@ end subroutine sampleCollision
   !!
   !! Perform implicit treatment
   !!
-  subroutine implicit(self, p, collDat, thisCycle, nextCycle)
+  subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)   :: self
     class(particle), intent(inout)       :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)   :: collDat
     class(particleDungeon),intent(inout) :: thisCycle
     class(particleDungeon),intent(inout) :: nextCycle
@@ -302,7 +312,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
       ! Store new sites in the next cycle dungeon
       r   = p % rGlobal()
 
-      do i=1,n
+      do i = 1,n
         call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG)
         dir = rotateVector(p % dirGlobal(), mu, phi)
 
@@ -321,6 +331,9 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
         call nextCycle % detain(pTemp)
         if (self % uniFissSites) call self % ufsField % storeFS(pTemp)
 
+        ! Report birth of new particle
+        call tally % reportSpawn(N_FISSION, p, pTemp)
+
       end do
     end if
 
@@ -346,9 +359,10 @@ end subroutine implicit
   !!
   !! Process capture reaction
   !!
-  subroutine capture(self, p, collDat, thisCycle, nextCycle)
+  subroutine capture(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)   :: self
     class(particle), intent(inout)       :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)   :: collDat
     class(particleDungeon),intent(inout) :: thisCycle
     class(particleDungeon),intent(inout) :: nextCycle
@@ -360,9 +374,10 @@ end subroutine capture
   !!
   !! Process fission reaction
   !!
-  subroutine fission(self, p, collDat, thisCycle, nextCycle)
+  subroutine fission(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)   :: self
     class(particle), intent(inout)       :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)   :: collDat
     class(particleDungeon),intent(inout) :: thisCycle
     class(particleDungeon),intent(inout) :: nextCycle
@@ -426,6 +441,9 @@ subroutine fission(self, p, collDat, thisCycle, nextCycle)
         call nextCycle % detain(pTemp)
         if (self % uniFissSites) call self % ufsField % storeFS(pTemp)
 
+        ! Report birth of new particle
+        call tally % reportSpawn(N_FISSION, p, pTemp)
+
       end do
     end if
 
@@ -438,9 +456,10 @@ end subroutine fission
   !!
   !! All CE elastic scattering happens in the CM frame
   !!
-  subroutine elastic(self, p, collDat, thisCycle, nextCycle)
+  subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)     :: self
     class(particle), intent(inout)         :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)     :: collDat
     class(particleDungeon),intent(inout)   :: thisCycle
     class(particleDungeon),intent(inout)   :: nextCycle
@@ -476,9 +495,10 @@ end subroutine elastic
   !!
   !! Process inelastic scattering
   !!
-  subroutine inelastic(self, p, collDat, thisCycle, nextCycle)
+  subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)     :: self
     class(particle), intent(inout)         :: p
+    type(tallyAdmin), intent(inout)        :: tally
     type(collisionData), intent(inout)     :: collDat
     class(particleDungeon),intent(inout)   :: thisCycle
     class(particleDungeon),intent(inout)   :: nextCycle
@@ -506,17 +526,19 @@ end subroutine inelastic
   !!
   !! Apply cutoffs
   !!
-  subroutine cutoffs(self, p, collDat, thisCycle, nextCycle)
+  subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle)
     class(neutronCEimp), intent(inout)   :: self
     class(particle), intent(inout)       :: p
+    type(tallyAdmin), intent(inout)      :: tally
     type(collisionData), intent(inout)   :: collDat
     class(particleDungeon),intent(inout) :: thisCycle
     class(particleDungeon),intent(inout) :: nextCycle
     real(defReal), dimension(3)          :: val
     real(defReal)                        :: minWgt, maxWgt, avWgt
 
     if (p % isDead) then
-    ! Do nothing
+
+      ! Do nothing !
 
     elseif (p % E < self % minE) then
       p % isDead = .true.
@@ -530,15 +552,18 @@ subroutine cutoffs(self, p, 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
-        call self % split(p, thisCycle, maxWgt)
+
+      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
     elseif ((self % splitting) .and. (p % w > self % maxWgt)) then
-      call self % split(p, thisCycle, self % maxWgt)
+      call self % split(p, tally, thisCycle, self % maxWgt)
 
     ! Roulette with fixed threshold and survival weight
     elseif ((self % roulette) .and. (p % w < self % minWgt)) then
@@ -567,24 +592,41 @@ end subroutine russianRoulette
   !!
   !! Split particle which has too large a weight
   !!
-  subroutine split(self, p, thisCycle, maxWgt)
+  subroutine split(self, p, tally, thisCycle, maxWgt)
     class(neutronCEimp), intent(inout)    :: self
     class(particle), intent(inout)        :: p
+    type(tallyAdmin), intent(inout)       :: tally
     class(particleDungeon), intent(inout) :: thisCycle
     real(defReal), intent(in)             :: maxWgt
+    type(particleState)                   :: pTemp
     integer(shortInt)                     :: mult, i
 
     ! This value must be at least 2
     mult = ceiling(p % w/maxWgt)
 
-    ! Decrease weight
-    p % w = p % w/mult
+    ! 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
 
     ! Add split particle's to the dungeon
-    do i = 1,mult-1
-      call thisCycle % detain(p)
+    do i = 1, mult-1
+      call thisCycle % detain(pTemp)
+      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
+
   end subroutine split
 
   !!