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 2 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 @@ -118,9 +118,11 @@ 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)
Expand Down
24 changes: 18 additions & 6 deletions CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ module neutronCEimp_class

use numPrecision
use endfConstants
use universalVariables, only : nameUFS, nameWW
use universalVariables, only : nameUFS, nameWW, maxSplit
use genericProcedures, only : fatalError, rotateVector, numToChar
use dictionary_class, only : dictionary
use RNG_class, only : RNG
Expand Down Expand Up @@ -190,8 +190,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 @@ -516,7 +516,8 @@ subroutine cutoffs(self, p, 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 @@ -530,7 +531,7 @@ 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
if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < maxSplit)) then
call self % split(p, thisCycle, maxWgt)
elseif (p % w < minWgt) then
call self % russianRoulette(p, avWgt)
Expand Down Expand Up @@ -572,19 +573,30 @@ subroutine split(self, p, thisCycle, maxWgt)
class(particle), intent(inout) :: p
class(particleDungeon), intent(inout) :: thisCycle
real(defReal), intent(in) :: maxWgt
integer(shortInt) :: mult, i
integer(shortInt) :: mult, i, splitCount

! This value must be at least 2
mult = ceiling(p % w/maxWgt)

! Limit maximum split
if (mult > maxSplit - p % splitCount + 1) then
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (mult > maxSplit - p % splitCount + 1) then
if (mult + p % splitCount > maxSplit) then

Maybe this condition will be a bit more readable?
I had a bit of trouble to work out the intention behind this if statement. It does limit number of 2ndary particles created by splitting so they do not exceed maxSplit, right? Or did I misunderstood something?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, sure! That was the intention, I might have made it more convoluted than it needed to be! Thanks

mult = maxSplit - p % splitCount + 1
end if

! Decrease weight
p % w = p % w/mult

! Save current particle splitCount
splitCount = p % splitCount

! Add split particle's to the dungeon
do i = 1,mult-1
p % splitCount = 0
call thisCycle % detain(p)
end do

p % splitCount = splitCount + mult

end subroutine split

!!
Expand Down
Loading
Loading