Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding MGimp processor and maximum number of split per particle #110

Merged
merged 6 commits into from
Jan 25, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions CollisionOperator/CollisionProcessors/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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))
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
38 changes: 28 additions & 10 deletions CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -89,6 +90,7 @@ module neutronCEimp_class
!! #impGen <logical>;#
!! #UFS <logical>;#
!! #weightWindows <logical>;#
!! #maxSplit <integer>;#
!! }
!!
type, public, extends(collisionProcessor) :: neutronCEimp
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down
Loading
Loading