Merge branch 'CambridgeNuclear:main' into main

valeriaRaffuzzi · Dec 12, 2023 · 5909f3f · 5909f3f
2 parents 2deb8cf + 7a5c199
commit 5909f3f
Show file tree

Hide file tree

Showing 35 changed files with 1,251 additions and 190 deletions.
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
@@ -1,10 +1,14 @@
 
 version: 2
 
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+
 sphinx:
   configuration: docs/conf.py
 
 python:
-  version: 3.7
   install:
-    - requirements: docs/requirements-rtd.txt
+    - requirements: docs/requirements-rtd.txt
diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90
@@ -120,7 +120,7 @@ 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)
+    call tally % reportInColl(p, .false.)
     call p % savePreCollision()
 
     ! Choose collision nuclide and general type (Scatter, Capture or Fission)
@@ -154,6 +154,9 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle)
     ! Apply post collision implicit treatments
     call self % cutoffs(p, collDat, thisCycle, nextCycle)
 
+    ! Update particle collision counter
+    p % collisionN = p % collisionN + 1
+
     ! Report out-of-collision
     call tally % reportOutColl(p, collDat % MT, collDat % muL)
 

diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
@@ -37,7 +37,7 @@ module neutronCEimp_class
 
   ! Scattering procedures
   use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, &
-                                     asymptoticInelasticScatter
+                                     asymptoticInelasticScatter, targetVelocity_DBRCXS
   implicit none
   private
 
@@ -63,6 +63,8 @@ module neutronCEimp_class
   !!              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)
+  !!  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)
   !!  roulette  -> roulettes particles below certain weight (off by defautl)
   !!  weightWindows -> uses a weight windows field (off by default)
@@ -88,9 +90,9 @@ module neutronCEimp_class
   type, public, extends(collisionProcessor) :: neutronCEimp
     private
     !! Nuclear Data block pointer -> public so it can be used by subclasses (protected member)
-    class(ceNeutronDatabase), pointer, public :: xsData => null()
-    class(ceNeutronMaterial), pointer, public :: mat    => null()
-    class(ceNeutronNuclide),  pointer, public :: nuc    => null()
+    class(ceNeutronDatabase), pointer, public :: xsData  => null()
+    class(ceNeutronMaterial), pointer, public :: mat     => null()
+    class(ceNeutronNuclide),  pointer, public :: nuc     => null()
     class(uniFissSitesField), pointer :: ufsField => null()
 
     !! Settings - private
@@ -101,6 +103,9 @@ module neutronCEimp_class
     real(defReal) :: avWgt
     real(defReal) :: thresh_E
     real(defReal) :: thresh_A
+    real(defReal) :: DBRCeMin
+    real(defReal) :: DBRCeMax
+
     ! Variance reduction options
     logical(defBool) :: weightWindows
     logical(defBool) :: splitting
@@ -175,6 +180,10 @@ subroutine init(self, dict)
     if( self % thresh_E < 0) call fatalError(Here,' -ve energyThreshold')
     if( self % thresh_A < 0) call fatalError(Here,' -ve massThreshold')
 
+    ! DBRC energy limits
+    call dict % getOrDefault(self % DBRCeMin,'DBRCeMin', (1.0E-8_defReal))
+    call dict % getOrDefault(self % DBRCeMax,'DBRCeMax', (200E-6_defReal))
+
     if (self % splitting) then
       if (self % maxWgt < 2 * self % minWgt) call fatalError(Here,&
               'Upper weight bound must be at least twice the lower weight bound')
@@ -231,7 +240,7 @@ subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
     collDat % nucIdx = self % mat % sampleNuclide(p % E, p % pRNG)
 
     self % nuc => ceNeutronNuclide_CptrCast(self % xsData % getNuclide(collDat % nucIdx))
-    if(.not.associated(self % mat)) call fatalError(Here, 'Failed to retive 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)
@@ -307,6 +316,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
         pTemp % dir = dir
         pTemp % E   = E_out
         pTemp % wgt = wgt
+        pTemp % collisionN = 0
 
         call nextCycle % detain(pTemp)
         if (self % uniFissSites) call self % ufsField % storeFS(pTemp)
@@ -411,6 +421,7 @@ subroutine fission(self, p, collDat, thisCycle, nextCycle)
         pTemp % dir = dir
         pTemp % E   = E_out
         pTemp % wgt = wgt
+        pTemp % collisionN = 0
 
         call nextCycle % detain(pTemp)
         if (self % uniFissSites) call self % ufsField % storeFS(pTemp)
@@ -434,7 +445,7 @@ subroutine elastic(self, p, collDat, thisCycle, nextCycle)
     class(particleDungeon),intent(inout)   :: thisCycle
     class(particleDungeon),intent(inout)   :: nextCycle
     class(uncorrelatedReactionCE), pointer :: reac
-    logical(defBool)                       :: isFixed
+    logical(defBool)                       :: isFixed, hasDBRC
     character(100),parameter :: Here = 'elastic (neutronCEimp_class.f90)'
 
     ! Get reaction
@@ -445,7 +456,11 @@ subroutine elastic(self, p, collDat, thisCycle, nextCycle)
     collDat % A =  self % nuc % getMass()
     collDat % kT = self % nuc % getkT()
 
-    isFixed = (p % E > collDat % kT * self % thresh_E) .and. (collDat % A > self % thresh_A)
+    ! 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)
 
     ! Apply criterion for Free-Gas vs Fixed Target scattering
     if (.not. reac % inCMFrame()) then
@@ -647,27 +662,54 @@ subroutine scatterFromMoving(self, p, collDat, reac)
     class(particle), intent(inout)             :: p
     type(collisionData),intent(inout)          :: collDat
     class(uncorrelatedReactionCE), intent(in)  :: reac
-    integer(shortInt)                          :: MT, nucIdx
+    class(ceNeutronNuclide), pointer           :: ceNuc0K
+    integer(shortInt)                          :: nucIdx
     real(defReal)                              :: A, kT, mu
     real(defReal),dimension(3)                 :: V_n           ! Neutron velocity (vector)
     real(defReal)                              :: U_n           ! Neutron speed (scalar)
     real(defReal),dimension(3)                 :: dir_pre       ! Pre-collision direction
     real(defReal),dimension(3)                 :: dir_post      ! Post-collicion direction
     real(defReal),dimension(3)                 :: V_t, V_cm     ! Target and CM velocity
     real(defReal)                              :: phi, dummy
+    real(defReal)                              :: maj
+    logical(defBool)                           :: inEnergyRange, hasDBRC
+    character(100), parameter :: Here = 'ScatterFromMoving (neutronCEimp_class.f90)'
 
-    ! Read data
-    MT     = collDat % MT
-    nucIdx = collDat % nucIdx
+    ! Read collision data
     A      = collDat % A
     kT     = collDat % kT
+    nucIdx = collDat % nucIdx
 
     ! Get neutron direction and velocity
     dir_pre = p % dirGlobal()
     V_n     = dir_pre * sqrt(p % E)
 
-    ! Sample velocity of target
-    V_t = targetVelocity_constXS(p % E, dir_pre, A, kT, p % pRNG)
+    ! Sample target velocity with constant XS or with DBRC
+    ! Check energy range
+    inEnergyRange = ((p % E <= self % DBRCeMax) .and. (self % DBRCeMin <= p % E))
+    ! Check if DBRC is on for this target nuclide
+    hasDBRC = self % nuc % hasDBRC()
+
+    if (inEnergyRange .and. hasDBRC) then
+
+      ! Retrieve 0K nuclide index from DBRC nuclide map
+      nucIdx = self % xsData % mapDBRCnuc % get(nucIdx)
+
+      ! 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')
+
+      ! Get elastic scattering 0K majorant
+      maj = self % xsData % getScattMicroMajXS(p % E, kT, A, nucIdx)
+
+      ! Use DBRC to sample target velocity
+      V_t = targetVelocity_DBRCXS(ceNuc0K, p % E, dir_pre, A, kT, p % pRNG, maj)
+
+    else
+      ! Constant cross section approximation
+      V_t = targetVelocity_constXS(p % E, dir_pre, A, kT, p % pRNG)
+
+    end if
 
     ! Calculate Centre-of-Mass velocity
     V_cm = (V_n + V_t *A)/(A+1)
@@ -694,6 +736,7 @@ subroutine scatterFromMoving(self, p, collDat, reac)
     p % E = U_n * U_n
     call p % point(dir_post)
     collDat % muL = dot_product(dir_pre, dir_post)
+
   end subroutine scatterFromMoving
 
 

diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90
@@ -31,7 +31,7 @@ module neutronCEstd_class
 
   ! Scattering procedures
   use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, &
-                                     asymptoticInelasticScatter
+                                     asymptoticInelasticScatter, targetVelocity_DBRCXS
   implicit none
   private
 
@@ -51,6 +51,8 @@ module neutronCEstd_class
   !!              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)
+  !!  DBRCeMin -> Minimum energy to which DBRC is applied
+  !!  DBRCeMax -> Maximum energy to which DBRC is applied
   !!
   !! Sample dictionary input:
   !!   collProcName {
@@ -73,6 +75,8 @@ module neutronCEstd_class
     real(defReal) :: maxE
     real(defReal) :: thresh_E
     real(defReal) :: thresh_A
+    real(defReal) :: DBRCeMin
+    real(defReal) :: DBRCeMax
 
   contains
     ! Initialisation procedure
@@ -122,6 +126,10 @@ subroutine init(self, dict)
     if( self % thresh_E < 0) call fatalError(Here,' -ve energyThreshold')
     if( self % thresh_A < 0) call fatalError(Here,' -ve massThreshold')
 
+    ! DBRC energy limits
+    call dict % getOrDefault(self % DBRCeMin,'DBRCeMin', (1.0E-8_defReal))
+    call dict % getOrDefault(self % DBRCeMax,'DBRCeMax', (200E-6_defReal))
+
   end subroutine init
 
   !!
@@ -154,7 +162,7 @@ subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
     collDat % nucIdx = self % mat % sampleNuclide(p % E, p % pRNG)
 
     self % nuc => ceNeutronNuclide_CptrCast(self % xsData % getNuclide(collDat % nucIdx))
-    if(.not.associated(self % mat)) call fatalError(Here, 'Failed to retive 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)
@@ -222,6 +230,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
         pTemp % dir = dir
         pTemp % E   = E_out
         pTemp % wgt = wgt
+        pTemp % collisionN = 0
 
         call nextCycle % detain(pTemp)
       end do
@@ -263,13 +272,13 @@ end subroutine fission
   !! All CE elastic scattering happens in the CM frame
   !!
   subroutine elastic(self, p, collDat, thisCycle, nextCycle)
-    class(neutronCEstd), intent(inout)   :: self
-    class(particle), intent(inout)       :: p
-    type(collisionData), intent(inout)   :: collDat
-    class(particleDungeon),intent(inout) :: thisCycle
-    class(particleDungeon),intent(inout) :: nextCycle
+    class(neutronCEstd), intent(inout)     :: self
+    class(particle), intent(inout)         :: p
+    type(collisionData), intent(inout)     :: collDat
+    class(particleDungeon),intent(inout)   :: thisCycle
+    class(particleDungeon),intent(inout)   :: nextCycle
     class(uncorrelatedReactionCE), pointer :: reac
-    logical(defBool)                       :: isFixed
+    logical(defBool)                       :: isFixed, hasDBRC
     character(100),parameter :: Here = 'elastic (neutronCEstd_class.f90)'
 
     ! Get reaction
@@ -280,7 +289,11 @@ subroutine elastic(self, p, collDat, thisCycle, nextCycle)
     collDat % A =  self % nuc % getMass()
     collDat % kT = self % nuc % getkT()
 
-    isFixed = (p % E > collDat % kT * self % thresh_E) .and. (collDat % A > self % thresh_A)
+    ! 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)
 
     ! Apply criterion for Free-Gas vs Fixed Target scattering
     if (.not. reac % inCMFrame()) then
@@ -404,24 +417,54 @@ subroutine scatterFromMoving(self, p, collDat, reac)
     class(particle), intent(inout)             :: p
     type(collisionData),intent(inout)          :: collDat
     class(uncorrelatedReactionCE), intent(in)  :: reac
+    class(ceNeutronNuclide), pointer           :: ceNuc0K
+    integer(shortInt)                          :: nucIdx
     real(defReal)                              :: A, kT, mu
     real(defReal),dimension(3)                 :: V_n           ! Neutron velocity (vector)
     real(defReal)                              :: U_n           ! Neutron speed (scalar)
     real(defReal),dimension(3)                 :: dir_pre       ! Pre-collision direction
     real(defReal),dimension(3)                 :: dir_post      ! Post-collicion direction
     real(defReal),dimension(3)                 :: V_t, V_cm     ! Target and CM velocity
     real(defReal)                              :: phi, dummy
+    real(defReal)                              :: maj
+    logical(defBool)                           :: inEnergyRange, hasDBRC
+    character(100), parameter :: Here = 'ScatterFromMoving (neutronCEstd_class.f90)'
 
-    ! Read data
+    ! Read collision data
     A      = collDat % A
     kT     = collDat % kT
+    nucIdx = collDat % nucIdx
 
     ! Get neutron direction and velocity
     dir_pre = p % dirGlobal()
     V_n     = dir_pre * sqrt(p % E)
 
-    ! Sample velocity of target
-    V_t = targetVelocity_constXS(p % E, dir_pre, A, kT, p % pRNG)
+    ! Sample target velocity with constant XS or with DBRC
+    ! Check energy range
+    inEnergyRange = ((p % E <= self % DBRCeMax) .and. (self % DBRCeMin <= p % E))
+    ! Check if DBRC is on for this target nuclide
+    hasDBRC = self % nuc % hasDBRC()
+
+    if (inEnergyRange .and. hasDBRC) then
+
+      ! Retrieve 0K nuclide index from DBRC nuclide map
+      nucIdx = self % xsData % mapDBRCnuc % get(nucIdx)
+
+      ! 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')
+
+      ! Get elastic scattering 0K majorant
+      maj = self % xsData % getScattMicroMajXS(p % E, kT, A, nucIdx)
+
+      ! Use DBRC to sample target velocity
+      V_t = targetVelocity_DBRCXS(ceNuc0K, p % E, dir_pre, A, kT, p % pRNG, maj)
+
+    else
+      ! Constant cross section approximation
+      V_t = targetVelocity_constXS(p % E, dir_pre, A, kT, p % pRNG)
+
+    end if
 
     ! Calculate Centre-of-Mass velocity
     V_cm = (V_n + V_t *A)/(A+1)

diff --git a/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90
@@ -175,6 +175,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
         pTemp % dir = dir
         pTemp % G   = G_out
         pTemp % wgt = wgt
+        pTemp % collisionN = 0
 
         call nextCycle % detain(pTemp)
       end do

diff --git a/CollisionOperator/collisionOperator_class.f90 b/CollisionOperator/collisionOperator_class.f90
@@ -114,7 +114,7 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
     class(particleDungeon),intent(inout)     :: thisCycle
     class(particleDungeon),intent(inout)     :: nextCycle
     integer(shortInt)                        :: idx, procType
-    character(100), parameter :: Here = 'collide ( collisionOperator_class.f90)'
+    character(100), parameter :: Here = 'collide (collisionOperator_class.f90)'
 
     ! Select processing index with ternary expression
     if(p % isMG) then