TMS

CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90

@@ -26,6 +26,7 @@ module collisionProcessor_inter
     real(defReal)     :: muL     = ONE  !! Cosine of deflection angle in LAB-frame
     real(defReal)     :: A       = ZERO !! Target Mass [Neutron Mass]
     real(defReal)     :: kT      = ZERO !! Target temperature [MeV]
+    real(defReal)     :: E       = ZERO !! Collision energy (could be relative to target) [MeV]
   end type
 
 
@@ -113,25 +114,35 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
     class(particleDungeon),intent(inout)     :: thisCycle
     class(particleDungeon),intent(inout)     :: nextCycle
     type(collisionData)                      :: collDat
-    character(100),parameter                 :: Here = ' collide (collisionProcessor.f90)'
+    logical(defBool)                         :: virtual
+    integer(shortInt)                        :: addCollision
+    character(100),parameter                 :: Here = 'collide (collisionProcessor.f90)'
 
     ! Load material index into data package
     collDat % matIdx = p % matIdx()
 
+    ! Choose collision nuclide and general type (Scatter, Capture or Fission)
+    call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)
+
+    ! In case of a TMS rejection, set collision as virtual
+    if (collDat % MT == noInteraction) then
+      virtual = .true.
+      addCollision = 0
+    else
+      virtual = .false.
+      addCollision = 1
+    end if
+
     ! 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
     ! tallies (e.g., collisionProbability) to return dubious results
-
-    call tally % reportInColl(p, .false.)
+    call tally % reportInColl(p, virtual)
 
     call p % savePreCollision()
 
-    ! Choose collision nuclide and general type (Scatter, Capture or Fission)
-    call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)
-
     ! Perform implicit treatment
-    call self % implicit(p, tally, collDat, thisCycle, nextCycle)
+    if (collDat % MT /= noInteraction) call self % implicit(p, tally, collDat, thisCycle, nextCycle)
 
     ! Select physics to be processed based on MT number
     select case(collDat % MT)
@@ -159,13 +170,13 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
     call self % cutoffs(p, tally, collDat, thisCycle, nextCycle)
 
     ! Update particle collision counter
-    p % collisionN = p % collisionN + 1
+    p % collisionN = p % collisionN + addCollision
 
     ! Report out-of-collision
     call tally % reportOutColl(p, collDat % MT, collDat % muL)
 
     ! Report end-of-history if particle was killed
-    if( p % isDead) then
+    if (p % isDead) then
       p % fate = ABS_FATE
       call tally % reportHist(p)
     end if

CollisionOperator/CollisionProcessors/neutronCEimp_class.f90

@@ -2,7 +2,7 @@ module neutronCEimp_class
 
   use numPrecision
   use endfConstants
-  use universalVariables,            only : nameUFS, nameWW
+  use universalVariables,            only : nameUFS, nameWW, REJECTED
   use genericProcedures,             only : fatalError, rotateVector, numToChar
   use dictionary_class,              only : dictionary
   use RNG_class,                     only : RNG
@@ -63,11 +63,11 @@ module neutronCEimp_class
   !!  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
+  !!  threshE  -> Energy threshold for explicit treatment of target nuclide movement [-].
+  !!              Target movement is sampled if neutron energy E < kT * threshE where
   !!              kT is target material temperature in [MeV]. (default = 400.0)
-  !!  thresh_A -> Mass threshold for explicit tratment of target nuclide movement [Mn].
-  !!              Target movment is sampled if target mass A < thresh_A. (default = 1.0)
+  !!  threshA  -> Mass threshold for explicit treatment of target nuclide movement [Mn].
+  !!              Target movment is sampled if target mass A < threshA. (default = 1.0)
   !!  DBRCeMin -> Minimum energy to which DBRC is applied
   !!  DBRCeMax -> Maximum energy to which DBRC is applied
   !!  splitting -> splits particles above certain weight (on by default)
@@ -107,8 +107,8 @@ module neutronCEimp_class
     real(defReal) :: minWgt
     real(defReal) :: maxWgt
     real(defReal) :: avWgt
-    real(defReal) :: thresh_E
-    real(defReal) :: thresh_A
+    real(defReal) :: threshE
+    real(defReal) :: threshA
     real(defReal) :: DBRCeMin
     real(defReal) :: DBRCeMax
     integer(shortInt) :: maxSplit
@@ -121,6 +121,7 @@ module neutronCEimp_class
     logical(defBool)  :: implicitSites ! Generates fission sites on every fissile collision
     logical(defBool)  :: uniFissSites
 
+    ! Variance reduction requirements
     type(weightWindowsField), pointer :: weightWindowsMap
 
   contains
@@ -166,8 +167,8 @@ subroutine init(self, dict)
     call dict % getOrDefault(self % maxE,'maxEnergy',20.0_defReal)
 
     ! Thermal scattering kernel thresholds
-    call dict % getOrDefault(self % thresh_E, 'energyThreshold', 400.0_defReal)
-    call dict % getOrDefault(self % thresh_A, 'massThreshold', 1.0_defReal)
+    call dict % getOrDefault(self % threshE, 'energyThreshold', 400.0_defReal)
+    call dict % getOrDefault(self % threshA, 'massThreshold', 1.0_defReal)
 
     ! Obtain settings for variance reduction
     call dict % getOrDefault(self % weightWindows,'weightWindows', .false.)
@@ -185,8 +186,8 @@ subroutine init(self, dict)
     if( self % minE < ZERO ) call fatalError(Here,'-ve minEnergy')
     if( self % maxE < ZERO ) call fatalError(Here,'-ve maxEnergy')
     if( self % minE >= self % maxE) call fatalError(Here,'minEnergy >= maxEnergy')
-    if( self % thresh_E < 0) call fatalError(Here,' -ve energyThreshold')
-    if( self % thresh_A < 0) call fatalError(Here,' -ve massThreshold')
+    if( self % threshE < 0) call fatalError(Here,' -ve energyThreshold')
+    if( self % threshA < 0) call fatalError(Here,' -ve massThreshold')
 
     ! DBRC energy limits
     call dict % getOrDefault(self % DBRCeMin,'DBRCeMin', (1.0E-8_defReal))
@@ -246,13 +247,19 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle)
     if(.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial')
 
     ! Select collision nuclide
-    collDat % nucIdx = self % mat % sampleNuclide(p % E, p % pRNG)
+    call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E)
+
+    ! If nuclide was rejected in TMS loop return to tracking
+    if (collDat % nucIdx == REJECTED) then
+      collDat % MT = noInteraction
+      return
+    end if
 
     self % nuc => ceNeutronNuclide_CptrCast(self % xsData % getNuclide(collDat % nucIdx))
-    if(.not.associated(self % mat)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide')
+    if (.not.associated(self % mat)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide')
 
     ! Select Main reaction channel
-    call self % nuc % getMicroXSs(microXss, p % E, p % pRNG)
+    call self % nuc % getMicroXSs(microXss, collDat % E, p % pRNG)
     r = p % pRNG % get()
     collDat % MT = microXss % invert(r)
 
@@ -281,13 +288,17 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)
 
     ! Generate fission sites if nuclide is fissile
     fiss_and_implicit = self % nuc % isFissile() .and. self % implicitSites
+
     if (fiss_and_implicit) then
+
       ! Obtain required data
       wgt   = p % w                ! Current weight
       k_eff = p % k_eff            ! k_eff for normalisation
       rand1 = p % pRNG % get()     ! Random number to sample sites
 
-      call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)
+      ! Retrieve cross section at the energy used for reaction sampling
+      call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)
+
       sig_nufiss = microXSs % nuFission
       sig_tot    = microXSs % total
 
@@ -339,19 +350,23 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)
 
     ! Perform implicit absorption
     if (self % implicitAbsorption) then
-      if(.not.fiss_and_implicit) then
-        call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)
+
+      if (.not.fiss_and_implicit) then
+        call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)
       end if
+
       sig_scatter  = microXSs % elasticScatter + microXSs % inelasticScatter
       sig_tot      = microXSs % total
       p % w        = p % w * sig_scatter/sig_tot
       ! Sample between elastic and inelastic
       totalElastic = microXSs % elasticScatter + microXSs % inelasticScatter
+
       if (p % pRNG % get() < microXSs % elasticScatter/totalElastic) then
         collDat % MT = N_N_elastic
       else
         collDat % MT = N_N_inelastic
       end if
+
     end if
 
   end subroutine implicit
@@ -391,12 +406,15 @@ subroutine fission(self, p, tally, collDat, thisCycle, nextCycle)
     character(100),parameter             :: Here = 'fission (neutronCEimp_class.f90)'
 
     if (.not.self % implicitSites) then
+
       ! Obtain required data
       wgt   = p % w                ! Current weight
       k_eff = p % k_eff            ! k_eff for normalisation
       rand1 = p % pRNG % get()     ! Random number to sample sites
 
-      call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)
+      ! Retrieve cross section at the energy used for reaction sampling
+      call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)
+
       sig_nufiss = microXSs % nuFission
       sig_fiss   = microXSs % fission
 
@@ -467,19 +485,28 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle)
     logical(defBool)                       :: isFixed, hasDBRC
     character(100),parameter :: Here = 'elastic (neutronCEimp_class.f90)'
 
+    ! Assess if thermal scattering data is needed or not
+    if (self % nuc % needsSabEl(p % E)) collDat % MT = N_N_ThermEL
+
     ! Get reaction
     reac => uncorrelatedReactionCE_CptrCast( self % xsData % getReaction(collDat % MT, collDat % nucIdx))
     if(.not.associated(reac)) call fatalError(Here,'Failed to get elastic neutron scatter')
 
     ! Scatter particle
     collDat % A =  self % nuc % getMass()
-    collDat % kT = self % nuc % getkT()
+
+    ! Retrieve kT from either material or nuclide
+    if (self % mat % useTMS(p % E)) then
+      collDat % kT = self % mat % kT
+    else
+      collDat % kT = self % nuc % getkT()
+    end if
 
     ! Check is DBRC is on
     hasDBRC = self % nuc % hasDBRC()
 
-    isFixed = (.not. hasDBRC) .and. (p % E > collDat % kT * self % thresh_E) &
-              & .and. (collDat % A > self % thresh_A)
+    isFixed = (.not. hasDBRC) .and. (p % E > collDat % kT * self % threshE) &
+              & .and. (collDat % A > self % threshA)
 
     ! Apply criterion for Free-Gas vs Fixed Target scattering
     if (.not. reac % inCMFrame()) then
@@ -506,7 +533,7 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle)
     character(100),parameter  :: Here =' inelastic (neutronCEimp_class.f90)'
 
     ! Invert inelastic scattering and Get reaction
-    collDat % MT = self % nuc % invertInelastic(p % E, p % pRNG)
+    collDat % MT = self % nuc % invertInelastic(collDat % E, p % pRNG)
     reac => uncorrelatedReactionCE_CptrCast( self % xsData % getReaction(collDat % MT, collDat % nucIdx))
     if(.not.associated(reac)) call fatalError(Here, "Failed to get scattering reaction")
 
@@ -680,7 +707,7 @@ subroutine scatterFromFixed(self, p, collDat, reac)
     ! Save incident energy
     E_out = p % E
 
-    if( MT == N_N_elastic) then
+    if (MT == N_N_elastic) then
       call asymptoticScatter(E_out, mu, collDat % A)
 
     else
@@ -739,7 +766,7 @@ subroutine scatterFromMoving(self, p, collDat, reac)
 
       ! Assign pointer for the 0K nuclide
       ceNuc0K => ceNeutronNuclide_CptrCast(self % xsData % getNuclide(nucIdx))
-      if(.not.associated(ceNuc0K)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide')
+      if (.not.associated(ceNuc0K)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide')
 
       ! Get elastic scattering 0K majorant
       maj = self % xsData % getScattMicroMajXS(p % E, kT, A, nucIdx)