Virtual collisions with surface tracking

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -59,7 +59,8 @@ module aceNeutronDatabase_class
   !!   nuclearData {
   !!   handles {
   !!   ce { type aceNeutronDatabase; DBRC (92238 94242); ures < 1 or 0 >;
-  !!        majorant < 1 or 0 >; aceLibrary <nuclear data path> ;} }
+  !!        majorant < 1 or 0 >; aceLibrary <nuclear data path> ;} 
+  !!        #avgDist 3.141;# }
   !!
   !! Public Members:
   !!   nuclides    -> array of aceNeutronNuclides with data
@@ -773,7 +774,8 @@ subroutine init(self, dict, ptr, silent )
     character(nameLen),dimension(:),allocatable      :: nucDBRC
     real(defReal)                                    :: A, nuckT, eUpSab, eUpSabNuc, &
                                                         eLowURR, eLowUrrNuc, alpha, &
-                                                        deltakT, eUpper, eLower, kT
+                                                        deltakT, eUpper, eLower, kT, &
+                                                        temp
     real(defReal), dimension(2)                      :: sabT
     integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0
     character(100), parameter :: Here = 'init (aceNeutronDatabase_class.f90)'
@@ -807,6 +809,18 @@ subroutine init(self, dict, ptr, silent )
       end do
     end do
 
+    ! Check for a minimum average collision distance
+    if (dict % isPresent('avgDist')) then
+      call dict % get(temp, 'avgDist')
+
+      if (temp <= ZERO) then
+        call fatalError(Here, 'Must have a finite, positive minimum average collision distance')
+      end if
+
+      self % collisionXS = ONE / temp
+
+    end if
+
     ! Get path to ACE library
     call dict % get(aceLibPath,'aceLibrary')
 

NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90

@@ -301,10 +301,14 @@ function getTrackingXS(self, p, matIdx, what) result(xs)
     select case(what)
 
       case (MATERIAL_XS)
-        xs = self % getTrackMatXS(p, matIdx)
+        if (matIdx == VOID_MAT) then
+          xs = self % collisionXS
+        else
+          xs = max(self % getTrackMatXS(p, matIdx), self % collisionXS)
+        end if
 
       case (MAJORANT_XS)
-        xs = self % getMajorantXS(p)
+        xs = max(self % getMajorantXS(p), self % collisionXS)
 
       case (TRACKING_XS)
 
@@ -355,6 +359,9 @@ function getTrackMatXS(self, p, matIdx) result(xs)
       print *,'Particle location: ', p % rGlobal()
       call fatalError(Here, 'Particle is in an undefined material with index: '&
               //numToChar(matIdx))
+    elseif (matIdx == VOID_MAT) then
+      xs = ZERO
+      return
     end if
 
     ! Check Cache and update if needed

NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90

@@ -45,6 +45,7 @@ module baseMgNeutronDatabase_class
   !!   nucData {
   !!     type baseMgNeutronDatabase;
   !!     PN P0;                        // or P1
+  !!     #avgDist 2.718; #
   !!   }
   !!
   !! Public Members:
@@ -103,10 +104,14 @@ function getTrackingXS(self, p, matIdx, what) result(xs)
     select case(what)
 
       case (MATERIAL_XS)
-        xs = self % getTrackMatXS(p, matIdx)
+        if (matIdx == VOID_MAT) then
+          xs = self % collisionXS
+        else
+          xs = max(self % getTrackMatXS(p, matIdx), self % collisionXS)
+        end if
 
       case (MAJORANT_XS)
-        xs = self % getMajorantXS(p)
+        xs = max(self % getMajorantXS(p), self % collisionXS)
 
       case (TRACKING_XS)
 
@@ -147,6 +152,9 @@ function getTrackMatXS(self, p, matIdx) result(xs)
       print *,'Particle location: ', p % rGlobal()
       call fatalError(Here, 'Particle is in an undefined material with index: '&
               //numToChar(matIdx))
+    else if (matIdx == VOID_MAT) then
+      xs = ZERO
+      return
     end if
 
     xs = self % getTotalMatXS(p, matIdx)
@@ -339,6 +347,7 @@ subroutine init(self, dict, ptr, silent)
     character(pathLen)                                 :: path
     character(nameLen)                                 :: scatterKey
     type(dictionary)                                   :: tempDict
+    real(defReal)                                      :: temp
     character(100), parameter :: Here = 'init (baseMgNeutronDatabase_class.f90)'
 
     ! Prevent reallocations
@@ -350,6 +359,18 @@ subroutine init(self, dict, ptr, silent)
     else
       loud = .true.
     end if
+
+    ! Check for a minimum average collision distance
+    if (dict % isPresent('avgDist')) then
+      call dict % get(temp, 'avgDist')
+
+      if (temp <= ZERO) then
+        call fatalError(Here, 'Must have a finite, positive minimum average collision distance')
+      end if
+
+      self % collisionXS = ONE / temp
+
+    end if
 
     ! Find number of materials and allocate space
     nMat = mm_nMat()

NuclearData/nuclearDatabase_inter.f90

@@ -20,6 +20,10 @@ module nuclearDatabase_inter
   !! So nuclear data objects for CE or MG neutron, photons or other particles are
   !! subclasses of this type.
   !!
+  !! Public Members:
+  !!   collisionXS -> XS determining the minimum average distance to collision (to be treated as virtual).
+  !!                  ZERO by default, i.e., flight distances can be infinite in void.
+  !!
   !! Interface:
   !!   getTrackingXS -> returns XS used to sample track length
   !!   getTrackMatXS -> returns material tracking xs, which could be different from the total (e.g., with TMS)
@@ -32,6 +36,7 @@ module nuclearDatabase_inter
   !!   kill          -> return to uninitialised state, clean memory
   !!
   type, public,abstract :: nuclearDatabase
+    real(defReal)                      :: collisionXS = ZERO
   contains
     procedure(init), deferred          :: init
     procedure(activate), deferred      :: activate

TransportOperator/transportOperatorHT_class.f90

@@ -156,24 +156,35 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle)
     class(particleDungeon),intent(inout)      :: thisCycle
     class(particleDungeon),intent(inout)      :: nextCycle
     integer(shortInt)                         :: event
-    real(defReal)                             :: sigmaT, dist
+    real(defReal)                             :: sigmaT, dist, sigmaTrack, invSigmaTrack
+    real(defReal), parameter                  :: tol  = 1.0E-12
     character(100), parameter :: Here = 'surfaceTracking (transportOperatorHT_class.f90)'
 
     STLoop: do
 
+      sigmaTrack = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS)
+
       ! Obtain the local cross-section, depending on the material
-      if (p % matIdx() == VOID_MAT) then
+      ! This branch is called in the case of voids with no imposed XS
+      if (sigmaTrack < tol) then
+
         dist = INFINITY
+        invSigmaTrack = INFINITY
+        sigmaT = ZERO
 
       else
-        sigmaT = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS)
-        dist = -log( p % pRNG % get()) / sigmaT
+
+        invSigmaTrack = ONE / sigmaTrack
+        dist = -log( p % pRNG % get()) * invSigmaTrack
+
+        ! Obtain the local cross-section
+        sigmaT = self % xsData % getTrackMatXS(p, p % matIdx())
 
         ! Should never happen! Catches NaN distances
         if (dist /= dist) call fatalError(Here, "Distance is NaN")
 
       end if
-
+      
       ! Save state before movement
       call p % savePrePath()
 
@@ -204,9 +215,18 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle)
           ! All is well
 
       end select
+
+      if (p % isDead) exit STLoop
 
-      ! Return if particle stoped at collision (not cell boundary)
-      if (event == COLL_EV .or. p % isDead) exit STLoop
+      ! Roll RNG to determine if the collision is real or virtual
+      ! Exit the loop if the collision is real, report collision if virtual
+      if (event == COLL_EV) then
+        if (p % pRNG % get() < sigmaT*invSigmaTrack) then
+          exit STLoop
+        else
+          call tally % reportInColl(p, .true.)
+        end if
+      end if
 
     end do STLoop
 

TransportOperator/transportOperatorST_class.f90

@@ -52,19 +52,30 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle)
     class(particleDungeon),intent(inout)      :: thisCycle
     class(particleDungeon),intent(inout)      :: nextCycle
     integer(shortInt)                         :: event
-    real(defReal)                             :: sigmaT, dist
+    real(defReal)                             :: sigmaT, dist, sigmaTrack, invSigmaTrack
     type(distCache)                           :: cache
+    real(defReal), parameter                  :: tol  = 1.0E-12
     character(100), parameter :: Here = 'surfaceTracking (transportOperatorST_class.f90)'
 
     STLoop: do
+
+      sigmaTrack = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS)
 
       ! Obtain the local cross-section, depending on the material
-      if (p % matIdx() == VOID_MAT) then
+      ! This branch is called in the case of voids with no imposed XS
+      if (sigmaTrack < tol) then
+
         dist = INFINITY
+        invSigmaTrack = INFINITY
+        sigmaT = ZERO
 
       else
-        sigmaT = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS)
-        dist = -log( p % pRNG % get()) / sigmaT
+
+        invSigmaTrack = ONE / sigmaTrack
+        dist = -log( p % pRNG % get()) * invSigmaTrack
+
+        ! Obtain the local cross-section
+        sigmaT = self % xsData % getTrackMatXS(p, p % matIdx())
 
         ! Should never happen! Catches NaN distances
         if (dist /= dist) call fatalError(Here, "Distance is NaN")
@@ -102,14 +113,23 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle)
         case(OVERLAP_MAT)
           print *, "Particle location: ", p % rGlobal()
           call fatalError(Here, "Particle is in overlapping cells")
-
+        
         case default
           ! All is well
 
       end select
 
-      ! Return if particle stopped at collision (not cell boundary)
-      if (event == COLL_EV .or. p % isDead) exit STLoop
+      if (p % isDead) exit STLoop
+
+      ! Roll RNG to determine if the collision is real or virtual
+      ! Exit the loop if the collision is real, report collision if virtual
+      if (event == COLL_EV) then
+        if (p % pRNG % get() < sigmaT*invSigmaTrack) then
+          exit STLoop
+        else
+          call tally % reportInColl(p, .true.)
+        end if
+      end if
 
     end do STLoop
 

docs/User Manual.rst

@@ -768,7 +768,10 @@ The **handles** definition is structured as the following: ::
       }
 
 The name of a handle has to be the same as defined in a ``physicsPackage`` under the
-keyword ``XSdata``.
+keyword ``XSdata``. The nuclear database can also be used to optionally set the minimum average
+collision distance for particles. This may be desirable in order to induce virtual collisions
+when using surface tracking in low density materials, for example. This can be done by using
+the ``avgDist`` keyword, followed by specifying the minimum average distance as desired.
 
 Otherwise, the possible **nuclear database** types allowed are:
 
@@ -786,11 +789,14 @@ from ACE files.
   to be applied.
 * majorant (*optional*, default = 1): 1 for true; 0 for false; flag to activate the
   pre-construction of a unionised majorant cross section
+* avgDist (*optional*, default = infinity): the minimum average distance until a
+  collision, which may be virtual. Used to obtain better statistics for the
+  collision estimator in low density materials, especially when using surface tracking.
 
 Example: ::
 
       ceData { type aceNuclearDatabase; aceLibrary ./myFolder/ACElib/JEF311.aceXS;
-      ures 1; DBRC (92238 94242)}
+      ures 1; DBRC (92238 94242); avgDist 32; }
 
 .. note::
    If DBRC is applied, the 0K cross section ace files of the relevant nuclides must
@@ -803,6 +809,9 @@ baseMgNeutronDatabase, used for multi-group data. In this case, the data is read
 from files provided by the user.
 
 * PN: includes a flag for anisotropy treatment. Could be ``P0`` or ``P1``
+* avgDist (*optional*, default = infinity): the minimum average distance until a
+  collision, which may be virtual. Used to obtain better statistics for the
+  collision estimator in low density materials, especially when using surface tracking.
 
 Example: ::