Skip to content

Commit

Permalink
Merge pull request #121 from valeriaRaffuzzi/TMS2
Browse files Browse the repository at this point in the history 
TMS
  • Loading branch information
@valeriaRaffuzzi
valeriaRaffuzzi authored Jul 26, 2024
2 parents 3dd3132 + 85aa16d commit 9078f8f
Show file tree
Hide file tree
Showing 52 changed files with 1,890 additions and 967 deletions.
29 changes: 20 additions & 9 deletions CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module collisionProcessor_inter
real(defReal) :: muL = ONE !! Cosine of deflection angle in LAB-frame
real(defReal) :: A = ZERO !! Target Mass [Neutron Mass]
real(defReal) :: kT = ZERO !! Target temperature [MeV]
real(defReal) :: E = ZERO !! Collision energy (could be relative to target) [MeV]
end type


Expand Down Expand Up @@ -113,25 +114,35 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
class(particleDungeon),intent(inout) :: thisCycle
class(particleDungeon),intent(inout) :: nextCycle
type(collisionData) :: collDat
character(100),parameter :: Here = ' collide (collisionProcessor.f90)'
logical(defBool) :: virtual
integer(shortInt) :: addCollision
character(100),parameter :: Here = 'collide (collisionProcessor.f90)'

! Load material index into data package
collDat % matIdx = p % matIdx()

! Choose collision nuclide and general type (Scatter, Capture or Fission)
call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)

! In case of a TMS rejection, set collision as virtual
if (collDat % MT == noInteraction) then
virtual = .true.
addCollision = 0
else
virtual = .false.
addCollision = 1
end if

! 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
! tallies (e.g., collisionProbability) to return dubious results

call tally % reportInColl(p, .false.)
call tally % reportInColl(p, virtual)

call p % savePreCollision()

! Choose collision nuclide and general type (Scatter, Capture or Fission)
call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle)

! Perform implicit treatment
call self % implicit(p, tally, collDat, thisCycle, nextCycle)
if (collDat % MT /= noInteraction) call self % implicit(p, tally, collDat, thisCycle, nextCycle)

! Select physics to be processed based on MT number
select case(collDat % MT)
Expand Down Expand Up @@ -159,13 +170,13 @@ subroutine collide(self, p, tally, thisCycle, nextCycle)
call self % cutoffs(p, tally, collDat, thisCycle, nextCycle)

! Update particle collision counter
p % collisionN = p % collisionN + 1
p % collisionN = p % collisionN + addCollision

! Report out-of-collision
call tally % reportOutColl(p, collDat % MT, collDat % muL)

! Report end-of-history if particle was killed
if( p % isDead) then
if (p % isDead) then
p % fate = ABS_FATE
call tally % reportHist(p)
end if
Expand Down
75 changes: 51 additions & 24 deletions CollisionOperator/CollisionProcessors/neutronCEimp_class.f90
View file
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, REJECTED
use genericProcedures, only : fatalError, rotateVector, numToChar
use dictionary_class, only : dictionary
use RNG_class, only : RNG
Expand Down Expand Up @@ -63,11 +63,11 @@ module neutronCEimp_class
!! 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
!! threshE -> Energy threshold for explicit treatment of target nuclide movement [-].
!! Target movement is sampled if neutron energy E < kT * threshE where
!! 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)
!! threshA -> Mass threshold for explicit treatment of target nuclide movement [Mn].
!! Target movment is sampled if target mass A < threshA. (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)
Expand Down Expand Up @@ -107,8 +107,8 @@ module neutronCEimp_class
real(defReal) :: minWgt
real(defReal) :: maxWgt
real(defReal) :: avWgt
real(defReal) :: thresh_E
real(defReal) :: thresh_A
real(defReal) :: threshE
real(defReal) :: threshA
real(defReal) :: DBRCeMin
real(defReal) :: DBRCeMax
integer(shortInt) :: maxSplit
Expand All @@ -121,6 +121,7 @@ module neutronCEimp_class
logical(defBool) :: implicitSites ! Generates fission sites on every fissile collision
logical(defBool) :: uniFissSites

! Variance reduction requirements
type(weightWindowsField), pointer :: weightWindowsMap

contains
Expand Down Expand Up @@ -166,8 +167,8 @@ subroutine init(self, dict)
call dict % getOrDefault(self % maxE,'maxEnergy',20.0_defReal)

! Thermal scattering kernel thresholds
call dict % getOrDefault(self % thresh_E, 'energyThreshold', 400.0_defReal)
call dict % getOrDefault(self % thresh_A, 'massThreshold', 1.0_defReal)
call dict % getOrDefault(self % threshE, 'energyThreshold', 400.0_defReal)
call dict % getOrDefault(self % threshA, 'massThreshold', 1.0_defReal)

! Obtain settings for variance reduction
call dict % getOrDefault(self % weightWindows,'weightWindows', .false.)
Expand All @@ -185,8 +186,8 @@ subroutine init(self, dict)
if( self % minE < ZERO ) call fatalError(Here,'-ve minEnergy')
if( self % maxE < ZERO ) call fatalError(Here,'-ve maxEnergy')
if( self % minE >= self % maxE) call fatalError(Here,'minEnergy >= maxEnergy')
if( self % thresh_E < 0) call fatalError(Here,' -ve energyThreshold')
if( self % thresh_A < 0) call fatalError(Here,' -ve massThreshold')
if( self % threshE < 0) call fatalError(Here,' -ve energyThreshold')
if( self % threshA < 0) call fatalError(Here,' -ve massThreshold')

! DBRC energy limits
call dict % getOrDefault(self % DBRCeMin,'DBRCeMin', (1.0E-8_defReal))
Expand Down Expand Up @@ -246,13 +247,19 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle)
if(.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial')

! Select collision nuclide
collDat % nucIdx = self % mat % sampleNuclide(p % E, p % pRNG)
call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E)

! If nuclide was rejected in TMS loop return to tracking
if (collDat % nucIdx == REJECTED) then
collDat % MT = noInteraction
return
end if

self % nuc => ceNeutronNuclide_CptrCast(self % xsData % getNuclide(collDat % nucIdx))
if(.not.associated(self % mat)) call fatalError(Here, 'Failed to retrieve 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)
call self % nuc % getMicroXSs(microXss, collDat % E, p % pRNG)
r = p % pRNG % get()
collDat % MT = microXss % invert(r)

Expand Down Expand Up @@ -281,13 +288,17 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)

! Generate fission sites if nuclide is fissile
fiss_and_implicit = self % nuc % isFissile() .and. self % implicitSites

if (fiss_and_implicit) then

! Obtain required data
wgt = p % w ! Current weight
k_eff = p % k_eff ! k_eff for normalisation
rand1 = p % pRNG % get() ! Random number to sample sites

call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)
! Retrieve cross section at the energy used for reaction sampling
call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)

sig_nufiss = microXSs % nuFission
sig_tot = microXSs % total

Expand Down Expand Up @@ -339,19 +350,23 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle)

! Perform implicit absorption
if (self % implicitAbsorption) then
if(.not.fiss_and_implicit) then
call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)

if (.not.fiss_and_implicit) then
call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)
end if

sig_scatter = microXSs % elasticScatter + microXSs % inelasticScatter
sig_tot = microXSs % total
p % w = p % w * sig_scatter/sig_tot
! Sample between elastic and inelastic
totalElastic = microXSs % elasticScatter + microXSs % inelasticScatter

if (p % pRNG % get() < microXSs % elasticScatter/totalElastic) then
collDat % MT = N_N_elastic
else
collDat % MT = N_N_inelastic
end if

end if

end subroutine implicit
Expand Down Expand Up @@ -391,12 +406,15 @@ subroutine fission(self, p, tally, collDat, thisCycle, nextCycle)
character(100),parameter :: Here = 'fission (neutronCEimp_class.f90)'

if (.not.self % implicitSites) then

! Obtain required data
wgt = p % w ! Current weight
k_eff = p % k_eff ! k_eff for normalisation
rand1 = p % pRNG % get() ! Random number to sample sites

call self % nuc % getMicroXSs(microXSs, p % E, p % pRNG)
! Retrieve cross section at the energy used for reaction sampling
call self % nuc % getMicroXSs(microXSs, collDat % E, p % pRNG)

sig_nufiss = microXSs % nuFission
sig_fiss = microXSs % fission

Expand Down Expand Up @@ -467,19 +485,28 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle)
logical(defBool) :: isFixed, hasDBRC
character(100),parameter :: Here = 'elastic (neutronCEimp_class.f90)'

! Assess if thermal scattering data is needed or not
if (self % nuc % needsSabEl(p % E)) collDat % MT = N_N_ThermEL

! Get reaction
reac => uncorrelatedReactionCE_CptrCast( self % xsData % getReaction(collDat % MT, collDat % nucIdx))
if(.not.associated(reac)) call fatalError(Here,'Failed to get elastic neutron scatter')

! Scatter particle
collDat % A = self % nuc % getMass()
collDat % kT = self % nuc % getkT()

! Retrieve kT from either material or nuclide
if (self % mat % useTMS(p % E)) then
collDat % kT = self % mat % kT
else
collDat % kT = self % nuc % getkT()
end if

! 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)
isFixed = (.not. hasDBRC) .and. (p % E > collDat % kT * self % threshE) &
& .and. (collDat % A > self % threshA)

! Apply criterion for Free-Gas vs Fixed Target scattering
if (.not. reac % inCMFrame()) then
Expand All @@ -506,7 +533,7 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle)
character(100),parameter :: Here =' inelastic (neutronCEimp_class.f90)'

! Invert inelastic scattering and Get reaction
collDat % MT = self % nuc % invertInelastic(p % E, p % pRNG)
collDat % MT = self % nuc % invertInelastic(collDat % E, p % pRNG)
reac => uncorrelatedReactionCE_CptrCast( self % xsData % getReaction(collDat % MT, collDat % nucIdx))
if(.not.associated(reac)) call fatalError(Here, "Failed to get scattering reaction")

Expand Down Expand Up @@ -680,7 +707,7 @@ subroutine scatterFromFixed(self, p, collDat, reac)
! Save incident energy
E_out = p % E

if( MT == N_N_elastic) then
if (MT == N_N_elastic) then
call asymptoticScatter(E_out, mu, collDat % A)

else
Expand Down Expand Up @@ -739,7 +766,7 @@ subroutine scatterFromMoving(self, p, collDat, reac)

! 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')
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)
Expand Down
Loading

0 comments on commit 9078f8f

Please sign in to comment.