From 13be97b0e5514652ef0de2588a084778544e2c51 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Thu, 30 Nov 2023 13:06:04 +0000 Subject: [PATCH 1/3] Adding MGimp processor and maximum number of split per particle --- .../CollisionProcessors/CMakeLists.txt | 5 +- .../collisionProcessorFactory_func.f90 | 8 +- .../collisionProcessor_inter.f90 | 3 +- .../neutronCEimp_class.f90 | 26 +- .../neutronMGimp_class.f90 | 371 ++++++++++++++++++ ParticleObjects/particle_class.f90 | 7 +- SharedModules/universalVariables.f90 | 5 +- 7 files changed, 412 insertions(+), 13 deletions(-) create mode 100644 CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 diff --git a/CollisionOperator/CollisionProcessors/CMakeLists.txt b/CollisionOperator/CollisionProcessors/CMakeLists.txt index ee923a0ed..8b8ab513a 100644 --- a/CollisionOperator/CollisionProcessors/CMakeLists.txt +++ b/CollisionOperator/CollisionProcessors/CMakeLists.txt @@ -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) diff --git a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 index 3f190f3a4..0ab4cb681 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 @@ -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 @@ -24,7 +25,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 @@ -59,6 +61,10 @@ subroutine new_collisionProcessor(new,dict) allocate(neutronMGstd :: new) call new % init(dict) + case('neutronMGimp') + allocate(neutronMGimp :: new) + call new % init(dict) + !*** NEW COLLISION PROCESSOR TEMPLATE ***! !case('') ! allocate( :: new) diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 index b19a106b7..a6fbc7014 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 @@ -118,8 +118,9 @@ 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) call p % savePreCollision() diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index 072e367bb..3e51c7378 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 @@ -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 @@ -181,8 +181,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 @@ -500,7 +500,10 @@ subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) real(defReal), dimension(3) :: val real(defReal) :: minWgt, maxWgt, avWgt - if (p % E < self % minE) then + if (p % isDead) then + ! Do nothing ! + + elseif (p % E < self % minE) then p % isDead = .true. ! Weight Windows treatment @@ -512,7 +515,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) @@ -554,19 +557,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 + 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 !! diff --git a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 new file mode 100644 index 000000000..3fa129cc8 --- /dev/null +++ b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 @@ -0,0 +1,371 @@ +module neutronMGimp_class + + use numPrecision + use endfConstants + use universalVariables, only : nameWW, maxSplit + use genericProcedures, only : fatalError, rotateVector, numToChar + use dictionary_class, only : dictionary + use RNG_class, only : RNG + + ! Particle types + use particle_class, only : particle, particleState, printType, P_NEUTRON + use particleDungeon_class, only : particleDungeon + + ! Abstract interface + use collisionProcessor_inter, only : collisionProcessor, collisionData ,init_super => init + + ! Nuclear Data Interface + use nuclearDataReg_mod, only : ndReg_getNeutronMG => getNeutronMG + use nuclearDatabase_inter, only : nuclearDatabase + use mgNeutronDatabase_inter, only : mgNeutronDatabase + use mgNeutronMaterial_inter, only : mgNeutronMaterial, mgNeutronMaterial_CptrCast + use reactionHandle_inter, only : reactionHandle + use multiScatterMG_class, only : multiScatterMG, multiScatterMG_CptrCast + use fissionMG_class, only : fissionMG, fissionMG_TptrCast + + ! Cross section packages + use neutronXsPackages_class, only : neutronMacroXSs + + ! Geometry and fields + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + use weightWindowsField_class, only : weightWindowsField, weightWindowsField_TptrCast + + + ! Nuclear Data + !use nuclearData_inter, only : nuclearData + !use perMaterialNuclearDataMG_inter, only : perMaterialNuclearDataMG + + ! Cross-section packages to interface with nuclear data + !use xsMacroSet_class, only : xsMacroSet, xsMacroSet_ptr + + implicit none + private + + !! + !! Standard (default) scalar collision processor for MG neutrons + !! -> Preforms implicit fission site generation + !! -> Preforms analog capture + !! -> Treats fission as capture (only implicit generation of 2nd-ary neutrons) + !! -> Does not create secondary non-neutron projectiles + !! + !! Settings: + !! NONE + !! + !! Sample dictionary input: + !! collProcName { + !! type neutronMGimp; + !! } + !! + type, public, extends(collisionProcessor) :: neutronMGimp + private + class(mgNeutronDatabase), pointer, public :: xsData => null() + class(mgNeutronMaterial), pointer, public :: mat => null() + + ! Variance reduction options + logical(defBool) :: weightWindows + type(weightWindowsField), pointer :: weightWindowsMap + + contains + ! Initialisation procedure + procedure :: init + + ! Implementation of customisable procedures + procedure :: sampleCollision + procedure :: implicit + procedure :: elastic + procedure :: inelastic + procedure :: capture + procedure :: fission + procedure :: cutoffs + + ! Variance reduction procedures + procedure, private :: split + procedure, private :: russianRoulette + + end type neutronMGimp + +contains + + !! + !! Initialise from dictionary + !! + subroutine init(self, dict) + class(neutronMGimp), intent(inout) :: self + class(dictionary), intent(in) :: dict + integer(shortInt) :: idx + character(100), parameter :: Here = 'init (neutronMGimp_class.f90)' + + ! Call superclass + call init_super(self, dict) + + ! Obtain settings for variance reduction + call dict % getOrDefault(self % weightWindows,'weightWindows', .false.) + + ! Sets up the weight windows field + if (self % weightWindows) then + idx = gr_fieldIdx(nameWW) + self % weightWindowsMap => weightWindowsField_TptrCast(gr_fieldPtr(idx)) + end if + + end subroutine init + + !! + !! Samples collision without any implicit treatment + !! + subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + type(neutronMacroXSs) :: macroXSs + real(defReal) :: r + character(100),parameter :: Here =' sampleCollision (neutronMGimp_class.f90)' + + ! Verify that particle is MG neutron + if( .not. p % isMG .or. p % type /= P_NEUTRON) then + call fatalError(Here, 'Supports only MG Neutron. Was given CE '//printType(p % type)) + end if + + ! Verify and load nuclear data pointer + self % xsData => ndReg_getNeutronMG() + if(.not.associated(self % xsData)) call fatalError(Here, "Failed to get active database for MG Neutron") + + ! Get and verify material pointer + self % mat => mgNeutronMaterial_CptrCast( self % xsData % getMaterial( p % matIdx())) + if(.not.associated(self % mat)) call fatalError(Here, "Failed to get MG Neutron Material") + + ! Select Main reaction channel + call self % mat % getMacroXSs(macroXSs, p % G, p % pRNG) + r = p % pRNG % get() + + collDat % MT = macroXSs % invert(r) + + end subroutine sampleCollision + + !! + !! Preform implicit treatment + !! + subroutine implicit(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + type(neutronMacroXSs) :: macroXSs + type(fissionMG),pointer :: fission + type(particleState) :: pTemp + real(defReal),dimension(3) :: r, dir + integer(shortInt) :: G_out, n, i + real(defReal) :: wgt, w0, rand1, mu, phi + real(defReal) :: sig_tot, k_eff, sig_nufiss + character(100),parameter :: Here = 'implicit (neutronMGimp_class.f90)' + + if ( self % mat % isFissile()) then + ! Obtain required data + wgt = p % w ! Current weight + w0 = p % preHistory % wgt ! Starting weight + k_eff = p % k_eff ! k_eff for normalisation + rand1 = p % pRNG % get() ! Random number to sample sites + + call self % mat % getMacroXSs(macroXSs, p % G, p % pRNG) + + sig_tot = macroXSs % total + sig_nuFiss = macroXSs % nuFission + + ! Sample number of fission sites generated + !n = int(wgt * sig_nuFiss/(sig_tot*k_eff) + r1, shortInt) + n = int(abs( (wgt * sig_nuFiss) / (w0 * sig_tot * k_eff)) + rand1, shortInt) + + ! Shortcut if no particles were samples + if (n < 1) return + + ! Get Fission reaction object + fission => fissionMG_TptrCast( self % xsData % getReaction(macroFission, collDat % matIdx)) + if (.not.associated(fission)) call fatalError(Here, 'Failed to getrive fissionMG reaction object') + + ! Store new sites in the next cycle dungeon + wgt = sign(w0, wgt) + r = p % rGlobal() + + do i=1,n + call fission % sampleOut(mu, phi, G_out, p % G, p % pRNG) + dir = rotateVector(p % dirGlobal(), mu, phi) + + ! Copy extra detail from parent particle (i.e. time, flags ect.) + pTemp = p + + ! Overwrite position, direction, energy group and weight + pTemp % r = r + pTemp % dir = dir + pTemp % G = G_out + pTemp % wgt = wgt + + call nextCycle % detain(pTemp) + end do + end if + + end subroutine implicit + + !! + !! Elastic Scattering + !! + subroutine elastic(self, p , collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + ! Do nothing. Should not be called + + end subroutine elastic + + !! + !! Preform scattering + !! + subroutine inelastic(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + class(multiScatterMG),pointer :: scatter + integer(shortInt) :: G_out ! Post-collision energy group + real(defReal) :: phi ! Azimuthal scatter angle + real(defReal) :: w_mul ! Weight multiplier + character(100),parameter :: Here = "inelastic (neutronMGimp_class.f90)" + + ! Assign MT number + collDat % MT = macroIEscatter + + ! Get Scatter object + scatter => multiScatterMG_CptrCast( self % xsData % getReaction(macroIEscatter, collDat % matIdx)) + if(.not.associated(scatter)) call fatalError(Here, "Failed to get scattering reaction object for MG neutron") + + ! Sample Mu and G_out + call scatter % sampleOut(collDat % muL, phi, G_out, p % G, p % pRNG) + + ! Read scattering multiplicity + w_mul = scatter % production(p % G, G_out) + + ! Update neutron state + p % G = G_out + p % w = p % w * w_mul + call p % rotate(collDat % muL, phi) + + end subroutine inelastic + + !! + !! Preform capture + !! + subroutine capture(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + p % isDead = .true. + + end subroutine capture + + !! + !! Preform fission + !! + subroutine fission(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + p % isDead = .true. + + end subroutine fission + + !! + !! Applay cutoffs or post-collision implicit treatment + !! + subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + real(defReal), dimension(3) :: val + real(defReal) :: minWgt, maxWgt, avWgt + + if (p % isDead) then + ! Do nothing ! + + ! Weight Windows treatment + elseif (self % weightWindows) then + val = self % weightWindowsMap % at(p) + minWgt = val(1) + maxWgt = val(2) + avWgt = val(3) + + ! 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) .and. (p % splitCount < maxSplit)) then + call self % split(p, thisCycle, maxWgt) + elseif (p % w < minWgt) then + call self % russianRoulette(p, avWgt) + end if + + end if + + end subroutine cutoffs + + !! + !! Perform Russian roulette on a particle + !! + subroutine russianRoulette(self, p, avWgt) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal), intent(in) :: avWgt + + if (p % pRNG % get() < (ONE - p % w/avWgt)) then + p % isDead = .true. + else + p % w = avWgt + end if + + end subroutine russianRoulette + + !! + !! Split particle which has too large a weight + !! + subroutine split(self, p, thisCycle, maxWgt) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + class(particleDungeon), intent(inout) :: thisCycle + real(defReal), intent(in) :: maxWgt + 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 + 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 + +end module neutronMGimp_class diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index 47556f487..2f1b34fe4 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -54,6 +54,7 @@ module particle_class integer(shortInt) :: matIdx = -1 ! Material index where particle is integer(shortInt) :: cellIdx = -1 ! Cell idx at the lowest coord level integer(shortInt) :: uniqueID = -1 ! Unique id at the lowest coord level + integer(shortInt) :: splitCount = 0 ! Counter of number of splits contains generic :: assignment(=) => fromParticle generic :: operator(.eq.) => equal_particleState @@ -112,7 +113,8 @@ module particle_class ! Particle processing information class(RNG), pointer :: pRNG => null() ! Pointer to RNG associated with the particle real(defReal) :: k_eff ! Value of default keff for implicit source generation - integer(shortInt) :: geomIdx ! Index of the geometry used by the particle + integer(shortInt) :: geomIdx ! Index of the geometry used by the particle + integer(shortInt) :: splitCount = 0 ! Counter of number of splits ! Archived snapshots of previous states type(particleState) :: preHistory @@ -269,6 +271,7 @@ pure subroutine particle_fromParticleState(LHS,RHS) LHS % isMG = RHS % isMG LHS % type = RHS % type LHS % time = RHS % time + LHS % splitCount = RHS % splitCount end subroutine particle_fromParticleState @@ -611,6 +614,7 @@ subroutine particleState_fromParticle(LHS,RHS) LHS % matIdx = RHS % coords % matIdx LHS % uniqueID = RHS % coords % uniqueId LHS % cellIdx = RHS % coords % cell() + LHS % splitCount = RHS % splitCount end subroutine particleState_fromParticle @@ -699,6 +703,7 @@ elemental subroutine kill_particleState(self) self % matIdx = -1 self % cellIdx = -1 self % uniqueID = -1 + self % splitCount = 0 end subroutine kill_particleState diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 9c133c1ee..81239996d 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -78,7 +78,8 @@ module universalVariables real(defReal), parameter :: joulesPerMeV = 1.60218e-13 ! Convert MeV to J ! Global name variables used to define specific geometry or field types - character(nameLen), parameter :: nameUFS = 'uniFissSites' - character(nameLen), parameter :: nameWW = 'WeightWindows' + character(nameLen), parameter :: nameUFS = 'uniFissSites' + character(nameLen), parameter :: nameWW = 'WeightWindows' + integer(shortInt), parameter :: maxSplit = 1000 end module universalVariables From 4fbf588d883477796c4b6a12c58f952767021755 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Thu, 25 Jan 2024 15:03:30 +0000 Subject: [PATCH 2/3] Making maxSplit a modifiable property --- .../neutronCEimp_class.f90 | 30 ++++++++++--------- .../neutronMGimp_class.f90 | 12 ++++---- SharedModules/universalVariables.f90 | 1 - 3 files changed, 23 insertions(+), 20 deletions(-) diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index 3e51c7378..e3dd68d5f 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 @@ -2,7 +2,7 @@ module neutronCEimp_class use numPrecision use endfConstants - use universalVariables, only : nameUFS, nameWW, maxSplit + use universalVariables, only : nameUFS, nameWW use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -102,12 +102,13 @@ module neutronCEimp_class real(defReal) :: thresh_E real(defReal) :: thresh_A ! 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 + integer(shortInt) :: maxSplit + 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 @@ -159,6 +160,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) @@ -515,7 +517,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) .and. (p % splitCount < maxSplit)) then + if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < self % maxSplit)) then call self % split(p, thisCycle, maxWgt) elseif (p % w < minWgt) then call self % russianRoulette(p, avWgt) @@ -563,23 +565,23 @@ subroutine split(self, p, thisCycle, maxWgt) mult = ceiling(p % w/maxWgt) ! Limit maximum split - if (mult > maxSplit - p % splitCount + 1) then - mult = maxSplit - p % splitCount + 1 + if (mult > self % maxSplit - p % splitCount + 1) then + mult = self % maxSplit - p % splitCount + 1 end if ! Decrease weight p % w = p % w/mult - + p % splitCount = p % splitCount + mult ! Save current particle splitCount - splitCount = p % splitCount + !splitCount = p % splitCount ! Add split particle's to the dungeon do i = 1,mult-1 - p % splitCount = 0 + !p % splitCount = 0 call thisCycle % detain(p) end do - p % splitCount = splitCount + mult + !p % splitCount = splitCount + mult end subroutine split diff --git a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 index 3fa129cc8..7560b4587 100644 --- a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 @@ -2,7 +2,7 @@ module neutronMGimp_class use numPrecision use endfConstants - use universalVariables, only : nameWW, maxSplit + use universalVariables, only : nameWW use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -62,7 +62,8 @@ module neutronMGimp_class class(mgNeutronMaterial), pointer, public :: mat => null() ! Variance reduction options - logical(defBool) :: weightWindows + integer(shortInt) :: maxSplit + logical(defBool) :: weightWindows type(weightWindowsField), pointer :: weightWindowsMap contains @@ -99,6 +100,7 @@ subroutine init(self, dict) call init_super(self, dict) ! Obtain settings for variance reduction + call dict % getOrDefault(self % maxSplit,'maxSplit', 1000) call dict % getOrDefault(self % weightWindows,'weightWindows', .false.) ! Sets up the weight windows field @@ -308,7 +310,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) .and. (p % splitCount < maxSplit)) then + if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < self % maxSplit)) then call self % split(p, thisCycle, maxWgt) elseif (p % w < minWgt) then call self % russianRoulette(p, avWgt) @@ -348,8 +350,8 @@ subroutine split(self, p, thisCycle, maxWgt) mult = ceiling(p % w/maxWgt) ! Limit maximum split - if (mult > maxSplit - p % splitCount + 1) then - mult = maxSplit - p % splitCount + 1 + if (mult > self % maxSplit - p % splitCount + 1) then + mult = self % maxSplit - p % splitCount + 1 end if ! Decrease weight diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 81239996d..bf3997446 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -80,6 +80,5 @@ module universalVariables ! Global name variables used to define specific geometry or field types character(nameLen), parameter :: nameUFS = 'uniFissSites' character(nameLen), parameter :: nameWW = 'WeightWindows' - integer(shortInt), parameter :: maxSplit = 1000 end module universalVariables From 3fc43ae7329b941da446b3634615cda8889e95b5 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Thu, 25 Jan 2024 16:16:44 +0000 Subject: [PATCH 3/3] Completing documentation --- .../neutronCEimp_class.f90 | 14 +++++++----- .../neutronMGimp_class.f90 | 22 +++++++++---------- .../neutronMGstd_class.f90 | 8 ------- ParticleObjects/particle_class.f90 | 4 ---- docs/User Manual.rst | 14 ++++++++++++ 5 files changed, 32 insertions(+), 30 deletions(-) diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index e3dd68d5f..98bd500ee 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 @@ -58,6 +58,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) @@ -83,6 +84,7 @@ module neutronCEimp_class !! #impGen ;# !! #UFS ;# !! #weightWindows ;# + !! #maxSplit ;# !! } !! type, public, extends(collisionProcessor) :: neutronCEimp @@ -101,8 +103,8 @@ module neutronCEimp_class real(defReal) :: avWgt real(defReal) :: thresh_E real(defReal) :: thresh_A - ! Variance reduction options integer(shortInt) :: maxSplit + ! Variance reduction options logical(defBool) :: weightWindows logical(defBool) :: splitting logical(defBool) :: roulette @@ -565,23 +567,23 @@ subroutine split(self, p, thisCycle, maxWgt) mult = ceiling(p % w/maxWgt) ! Limit maximum split - if (mult > self % maxSplit - p % splitCount + 1) then + if (mult + p % splitCount > self % maxSplit) then mult = self % maxSplit - p % splitCount + 1 end if ! Decrease weight p % w = p % w/mult - p % splitCount = p % splitCount + mult + ! Save current particle splitCount - !splitCount = p % splitCount + splitCount = p % splitCount ! Add split particle's to the dungeon do i = 1,mult-1 - !p % splitCount = 0 + p % splitCount = 0 call thisCycle % detain(p) end do - !p % splitCount = splitCount + mult + p % splitCount = splitCount + mult end subroutine split diff --git a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 index 7560b4587..b685b15c5 100644 --- a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 @@ -30,30 +30,26 @@ module neutronMGimp_class use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr use weightWindowsField_class, only : weightWindowsField, weightWindowsField_TptrCast - - ! Nuclear Data - !use nuclearData_inter, only : nuclearData - !use perMaterialNuclearDataMG_inter, only : perMaterialNuclearDataMG - - ! Cross-section packages to interface with nuclear data - !use xsMacroSet_class, only : xsMacroSet, xsMacroSet_ptr - implicit none private !! - !! Standard (default) scalar collision processor for MG neutrons + !! Scalar collision processor for MG neutrons !! -> Preforms implicit fission site generation !! -> Preforms analog capture !! -> Treats fission as capture (only implicit generation of 2nd-ary neutrons) !! -> Does not create secondary non-neutron projectiles + !! -> Supports the use of weight windows !! !! Settings: - !! NONE + !! weightWindows -> uses a weight windows field (off by default) + !! maxSplit -> maximum number of splits allowed per particle (default = 1000) !! !! Sample dictionary input: !! collProcName { !! type neutronMGimp; + !! #weightWindows ;# + !! #maxSplit ;# !! } !! type, public, extends(collisionProcessor) :: neutronMGimp @@ -61,8 +57,10 @@ module neutronMGimp_class class(mgNeutronDatabase), pointer, public :: xsData => null() class(mgNeutronMaterial), pointer, public :: mat => null() - ! Variance reduction options + !! Settings - private integer(shortInt) :: maxSplit + + ! Variance reduction options logical(defBool) :: weightWindows type(weightWindowsField), pointer :: weightWindowsMap @@ -350,7 +348,7 @@ subroutine split(self, p, thisCycle, maxWgt) mult = ceiling(p % w/maxWgt) ! Limit maximum split - if (mult > self % maxSplit - p % splitCount + 1) then + if (mult + p % splitCount > self % maxSplit) then mult = self % maxSplit - p % splitCount + 1 end if diff --git a/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 index b357319cb..019b694ce 100644 --- a/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 @@ -25,14 +25,6 @@ module neutronMGstd_class ! Cross section packages use neutronXsPackages_class, only : neutronMacroXSs - - ! Nuclear Data - !use nuclearData_inter, only : nuclearData - !use perMaterialNuclearDataMG_inter, only : perMaterialNuclearDataMG - - ! Cross-section packages to interface with nuclear data - !use xsMacroSet_class, only : xsMacroSet, xsMacroSet_ptr - implicit none private diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index 2f1b34fe4..3a295cc66 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -54,7 +54,6 @@ module particle_class integer(shortInt) :: matIdx = -1 ! Material index where particle is integer(shortInt) :: cellIdx = -1 ! Cell idx at the lowest coord level integer(shortInt) :: uniqueID = -1 ! Unique id at the lowest coord level - integer(shortInt) :: splitCount = 0 ! Counter of number of splits contains generic :: assignment(=) => fromParticle generic :: operator(.eq.) => equal_particleState @@ -271,7 +270,6 @@ pure subroutine particle_fromParticleState(LHS,RHS) LHS % isMG = RHS % isMG LHS % type = RHS % type LHS % time = RHS % time - LHS % splitCount = RHS % splitCount end subroutine particle_fromParticleState @@ -614,7 +612,6 @@ subroutine particleState_fromParticle(LHS,RHS) LHS % matIdx = RHS % coords % matIdx LHS % uniqueID = RHS % coords % uniqueId LHS % cellIdx = RHS % coords % cell() - LHS % splitCount = RHS % splitCount end subroutine particleState_fromParticle @@ -703,7 +700,6 @@ elemental subroutine kill_particleState(self) self % matIdx = -1 self % cellIdx = -1 self % uniqueID = -1 - self % splitCount = 0 end subroutine kill_particleState diff --git a/docs/User Manual.rst b/docs/User Manual.rst index e1512e501..3ea998af8 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -278,6 +278,7 @@ neutronCEimp, to perform implicit collision processing * minWgt (*optional*, default = 0.25): minimum particle weight for rouletting * maxWgt (*optional*, default = 1.25): maximum particle weight for splitting * avgWgt (*optional*, default = 0.5): weight of a particle on surviving rouletting +* maxSplit (*optional*, default = 1000): maximum number of splits allowed per particle * impAbs (*optional*, default = 0): 1 for true; 0 for false; enables implicit capture * impGen (*optional*, default = 1): 1 for true; 0 for false; enables implicit fission sites generation @@ -302,6 +303,19 @@ Example: :: collisionOperator { neutronMG { type neutronMGstd; } } +neutronMGimp +############ + +neutronMGimp, to perform implicit collision processing + +* maxSplit (*optional*, default = 1000): maximum number of splits allowed per particle +* weightWindows (*optional*, default = 0): 1 for true; 0 for false; enables the use of + weight windows + +Example: :: + + collisionOperator { neutronMG { type neutronMGimp; weightWindows 1; maxSplit 50; } } + Weight Windows --------------