Skip to content

Commit

Permalink
Merge pull request #108 from valeriaRaffuzzi/fixedSourceMgXss
Browse files Browse the repository at this point in the history 
Adding spawn report to tally Admin
  • Loading branch information
@valeriaRaffuzzi
valeriaRaffuzzi authored Jan 25, 2024
2 parents 6c2ba03 + cbc5785 commit 58df62d
Show file tree
Hide file tree
Showing 16 changed files with 282 additions and 1,607 deletions.
22 changes: 12 additions & 10 deletions CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,13 +86,15 @@ module collisionProcessor_inter
!! Procedure interface for all customisable actions associated with
!! processing of sollision event (scatter, fission etc.)
!!
subroutine collisionAction(self, p, collDat, thisCycle, nextCycle)
subroutine collisionAction(self, p, tally, collDat, thisCycle, nextCycle)
import :: collisionProcessor, &
collisionData, &
tallyAdmin, &
particle,&
particleDungeon
class(collisionProcessor), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand All @@ -104,7 +106,7 @@ end subroutine collisionAction
!!
!! Generic flow of collision processing
!!
subroutine collide(self, p, tally ,thisCycle, nextCycle)
subroutine collide(self, p, tally, thisCycle, nextCycle)
class(collisionProcessor), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
Expand All @@ -118,30 +120,30 @@ 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)
call self % sampleCollision(p, collDat, thisCycle, nextCycle)
call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)

! Perform implicit treatment
call self % implicit(p, collDat, thisCycle, nextCycle)
call self % implicit(p, tally, collDat, thisCycle, nextCycle)

! Select physics to be processed based on MT number
select case(collDat % MT)
case(N_N_elastic, macroAllScatter)
call self % elastic(p, collDat, thisCycle, nextCycle)
call self % elastic(p, tally, collDat, thisCycle, nextCycle)

case(N_N_inelastic, macroIEScatter)
call self % inelastic(p, collDat, thisCycle, nextCycle)
call self % inelastic(p, tally, collDat, thisCycle, nextCycle)

case(N_DISAP, macroCapture)
call self % capture(p, collDat, thisCycle, nextCycle)
call self % capture(p, tally, collDat, thisCycle, nextCycle)

case(N_FISSION, macroFission)
call self % fission(p, collDat, thisCycle, nextCycle)
call self % fission(p, tally, collDat, thisCycle, nextCycle)

case(noInteraction)
! Do nothing
Expand All @@ -152,7 +154,7 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle)
end select

! Apply post collision implicit treatments
call self % cutoffs(p, collDat, thisCycle, nextCycle)
call self % cutoffs(p, tally, collDat, thisCycle, nextCycle)

! Update particle collision counter
p % collisionN = p % collisionN + 1
Expand Down
54 changes: 39 additions & 15 deletions CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,10 @@ module neutronCEimp_class
! Scattering procedures
use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, &
asymptoticInelasticScatter, targetVelocity_DBRCXS

! Tally interfaces
use tallyAdmin_class, only : tallyAdmin

implicit none
private

Expand Down Expand Up @@ -213,9 +217,10 @@ end subroutine init
!!
!! Samples collision without any implicit treatment
!!
subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -252,9 +257,10 @@ end subroutine sampleCollision
!!
!! Perform implicit treatment
!!
subroutine implicit(self, p, collDat, thisCycle, nextCycle)
subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -302,7 +308,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
! Store new sites in the next cycle dungeon
r = p % rGlobal()

do i=1,n
do i = 1,n
call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG)
dir = rotateVector(p % dirGlobal(), mu, phi)

Expand All @@ -321,6 +327,9 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
call nextCycle % detain(pTemp)
if (self % uniFissSites) call self % ufsField % storeFS(pTemp)

! Report birth of new particle
call tally % reportSpawn(N_FISSION, p, pTemp)

end do
end if

Expand All @@ -346,9 +355,10 @@ end subroutine implicit
!!
!! Process capture reaction
!!
subroutine capture(self, p, collDat, thisCycle, nextCycle)
subroutine capture(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand All @@ -360,9 +370,10 @@ end subroutine capture
!!
!! Process fission reaction
!!
subroutine fission(self, p, collDat, thisCycle, nextCycle)
subroutine fission(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -426,6 +437,9 @@ subroutine fission(self, p, collDat, thisCycle, nextCycle)
call nextCycle % detain(pTemp)
if (self % uniFissSites) call self % ufsField % storeFS(pTemp)

! Report birth of new particle
call tally % reportSpawn(N_FISSION, p, pTemp)

end do
end if

Expand All @@ -438,9 +452,10 @@ end subroutine fission
!!
!! All CE elastic scattering happens in the CM frame
!!
subroutine elastic(self, p, collDat, thisCycle, nextCycle)
subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -476,9 +491,10 @@ end subroutine elastic
!!
!! Process inelastic scattering
!!
subroutine inelastic(self, p, collDat, thisCycle, nextCycle)
subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -506,9 +522,10 @@ end subroutine inelastic
!!
!! Apply cutoffs
!!
subroutine cutoffs(self, p, collDat, thisCycle, nextCycle)
subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand All @@ -531,14 +548,14 @@ 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
call self % split(p, thisCycle, maxWgt)
call self % split(p, tally, thisCycle, maxWgt)
elseif (p % w < minWgt) then
call self % russianRoulette(p, avWgt)
end if

! Splitting with fixed threshold
elseif ((self % splitting) .and. (p % w > self % maxWgt)) then
call self % split(p, thisCycle, self % maxWgt)
call self % split(p, tally, thisCycle, self % maxWgt)

! Roulette with fixed threshold and survival weight
elseif ((self % roulette) .and. (p % w < self % minWgt)) then
Expand Down Expand Up @@ -567,24 +584,31 @@ end subroutine russianRoulette
!!
!! Split particle which has too large a weight
!!
subroutine split(self, p, thisCycle, maxWgt)
subroutine split(self, p, tally, thisCycle, maxWgt)
class(neutronCEimp), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
class(particleDungeon), intent(inout) :: thisCycle
real(defReal), intent(in) :: maxWgt
type(particleState) :: pTemp
integer(shortInt) :: mult, i

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

! Decrease weight
p % w = p % w/mult
! Copy particle to a particle state
pTemp = p
pTemp % wgt = p % w/mult

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

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

end subroutine split

!!
Expand Down
31 changes: 23 additions & 8 deletions CollisionOperator/CollisionProcessors/neutronCEstd_class.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,10 @@ module neutronCEstd_class
! Scattering procedures
use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, &
asymptoticInelasticScatter, targetVelocity_DBRCXS

! Tally interfaces
use tallyAdmin_class, only : tallyAdmin

implicit none
private

Expand Down Expand Up @@ -135,9 +139,10 @@ end subroutine init
!!
!! Samples collision without any implicit treatment
!!
subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -174,9 +179,10 @@ end subroutine sampleCollision
!!
!! Perform implicit treatment
!!
subroutine implicit(self, p, collDat, thisCycle, nextCycle)
subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -216,7 +222,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
wgt = sign(w0, wgt)
r = p % rGlobal()

do i=1,n
do i = 1,n
call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG)
dir = rotateVector(p % dirGlobal(), mu, phi)

Expand All @@ -233,6 +239,10 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle)
pTemp % collisionN = 0

call nextCycle % detain(pTemp)

! Report birth of new particle
call tally % reportSpawn(N_FISSION, p, pTemp)

end do
end if

Expand All @@ -241,9 +251,10 @@ end subroutine implicit
!!
!! Process capture reaction
!!
subroutine capture(self, p, collDat, thisCycle, nextCycle)
subroutine capture(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand All @@ -255,9 +266,10 @@ end subroutine capture
!!
!! Process fission reaction
!!
subroutine fission(self, p, collDat, thisCycle, nextCycle)
subroutine fission(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand All @@ -271,9 +283,10 @@ end subroutine fission
!!
!! All CE elastic scattering happens in the CM frame
!!
subroutine elastic(self, p, collDat, thisCycle, nextCycle)
subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -309,9 +322,10 @@ end subroutine elastic
!!
!! Process inelastic scattering
!!
subroutine inelastic(self, p, collDat, thisCycle, nextCycle)
subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down Expand Up @@ -339,9 +353,10 @@ end subroutine inelastic
!!
!! Apply cutoffs
!!
subroutine cutoffs(self, p, collDat, thisCycle, nextCycle)
subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle)
class(neutronCEstd), intent(inout) :: self
class(particle), intent(inout) :: p
type(tallyAdmin), intent(inout) :: tally
type(collisionData), intent(inout) :: collDat
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
Expand Down
Loading

0 comments on commit 58df62d

Please sign in to comment.