From 85aa16da043137983210989db5628215ef2ce5ba Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi" Date: Fri, 26 Jul 2024 12:54:17 +0100 Subject: [PATCH] Implement TMS --- .../collisionProcessor_inter.f90 | 29 +- .../neutronCEimp_class.f90 | 75 ++- .../neutronCEstd_class.f90 | 91 +-- CollisionOperator/scatteringKernels_func.f90 | 86 ++- InputFiles/Benchmarks/BEAVRS/BEAVRS | 382 ++++++------ InputFiles/Benchmarks/BEAVRS/BEAVRS2D | 382 ++++++------ InputFiles/Benchmarks/HoogenboomMartin/HMcore | 12 - InputFiles/Benchmarks/MCNP/BIGTEN | 4 +- InputFiles/Benchmarks/MCNP/Falstaff | 6 +- InputFiles/Benchmarks/MCNP/IEUMF04 | 16 +- InputFiles/Benchmarks/MCNP/Jezebel | 2 +- InputFiles/Benchmarks/MCNP/Jezebel233 | 2 +- InputFiles/Benchmarks/MCNP/Jezebel240 | 2 +- InputFiles/Benchmarks/MCNP/LEUST02 | 4 +- InputFiles/Benchmarks/MCNP/ORNL10 | 4 +- InputFiles/Benchmarks/MCNP/ORNL11 | 4 +- InputFiles/Benchmarks/MCNP/PNL2 | 2 +- InputFiles/Benchmarks/MCNP/PUMF11 | 4 +- InputFiles/Benchmarks/MCNP/SB212 | 8 +- InputFiles/Benchmarks/MCNP/Tinkertoy2 | 10 +- InputFiles/Benchmarks/MCNP/U233MF05 | 4 +- .../Benchmarks/SourceConvergence/checkerboard | 7 +- .../Tests/aceNeutronDatabase_iTest.f90 | 7 +- .../Tests/thermalScatteringData_iTest.f90 | 7 +- .../Tests/urrProbabilityTables_iTest.f90 | 7 +- .../aceDatabase/aceNeutronDatabase_class.f90 | 564 +++++++++++++----- .../aceDatabase/aceNeutronNuclide_class.f90 | 129 +++- .../ceNeutronData/ceNeutronCache_mod.f90 | 47 +- .../ceNeutronData/ceNeutronDatabase_inter.f90 | 116 +++- .../ceNeutronData/ceNeutronMaterial_class.f90 | 350 ++++++++--- .../ceNeutronData/ceNeutronNuclide_inter.f90 | 44 +- NuclearData/materialMenu_mod.f90 | 21 +- .../baseMgNeutronDatabase_class.f90 | 24 +- NuclearData/nuclearDatabase_inter.f90 | 34 +- .../testNeutronDatabase_class.f90 | 18 +- ParticleObjects/particle_class.f90 | 7 +- SharedModules/numPrecision.f90 | 7 +- SharedModules/universalVariables.f90 | 8 +- .../TallyClerks/Tests/collisionClerk_test.f90 | 4 +- .../Tests/keffImplicitClerk_test.f90 | 1 + Tallies/TallyClerks/Tests/mgXsClerk_test.f90 | 2 +- Tallies/TallyClerks/collisionClerk_class.f90 | 47 +- .../collisionProbabilityClerk_class.f90 | 7 +- .../TallyClerks/keffImplicitClerk_class.f90 | 33 +- Tallies/TallyClerks/mgXsClerk_class.f90 | 92 +-- Tallies/TallyClerks/simpleFMClerk_class.f90 | 61 +- .../TallyResponses/macroResponse_class.f90 | 8 +- .../TallyResponses/microResponse_class.f90 | 9 +- .../transportOperatorDT_class.f90 | 3 +- .../transportOperatorHT_class.f90 | 7 +- TransportOperator/transportOperator_inter.f90 | 2 +- docs/User Manual.rst | 55 +- 52 files changed, 1890 insertions(+), 967 deletions(-) diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 index 466951106..f14c3da75 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 @@ -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 @@ -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) @@ -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 diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index 89d5a09d7..416d36a2c 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, REJECTED use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -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) @@ -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 @@ -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 @@ -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.) @@ -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)) @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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") @@ -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 @@ -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) diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 index e2ddf2f44..e3544b64e 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 @@ -2,6 +2,7 @@ module neutronCEstd_class use numPrecision use endfConstants + use universalVariables, only : REJECTED use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -31,7 +32,8 @@ module neutronCEstd_class ! Scattering procedures use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, & - asymptoticInelasticScatter, targetVelocity_DBRCXS + asymptoticInelasticScatter, targetVelocity_DBRCXS, & + relativeEnergy_constXS ! Tally interfaces use tallyAdmin_class, only : tallyAdmin @@ -50,11 +52,11 @@ module neutronCEstd_class !! minE -> minimum energy cut-off [MeV] (default = 1.0E-11) !! maxE -> maximum energy. Higher energies are set to maximum (not re-rolled) [MeV] !! (default = 20.0) - !! 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) - !! 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) + !! 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) + !! threshA -> Mass threshold for explicit treatment of target nuclide movement [Mn]. + !! Target movement 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 !! @@ -77,8 +79,8 @@ module neutronCEstd_class !! Settings - private real(defReal) :: minE real(defReal) :: maxE - real(defReal) :: thresh_E - real(defReal) :: thresh_A + real(defReal) :: threshE + real(defReal) :: threshA real(defReal) :: DBRCeMin real(defReal) :: DBRCeMax @@ -99,6 +101,7 @@ module neutronCEstd_class procedure,private :: scatterFromFixed procedure,private :: scatterFromMoving procedure,private :: scatterInLAB + end type neutronCEstd contains @@ -120,15 +123,15 @@ 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) ! Verify settings - 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 % 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 % 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)) @@ -151,26 +154,32 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) character(100),parameter :: Here = 'sampleCollision (neutronCEstd_class.f90)' ! Verify that particle is CE neutron - if(p % isMG .or. p % type /= P_NEUTRON) then + if (p % isMG .or. p % type /= P_NEUTRON) then call fatalError(Here, 'Supports only CE Neutron. Was given MG '//printType(p % type)) end if ! Verify and load nuclear data pointer self % xsData => ndReg_getNeutronCE() - if(.not.associated(self % xsData)) call fatalError(Here, 'There is no active Neutron CE data!') + if (.not.associated(self % xsData)) call fatalError(Here, 'There is no active Neutron CE data!') ! Verify and load material pointer - self % mat => ceNeutronMaterial_CptrCast( self % xsData % getMaterial( p % matIdx())) - if(.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial') + self % mat => ceNeutronMaterial_CptrCast(self % xsData % getMaterial(p % matIdx())) + 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 % nuc)) 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) @@ -196,14 +205,17 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) character(100),parameter :: Here = 'implicit (neutronCEstd_class.f90)' ! Generate fission sites if nuclide is fissile - if ( self % nuc % isFissile()) then + if (self % nuc % 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 % 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 @@ -294,19 +306,28 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) logical(defBool) :: isFixed, hasDBRC character(100),parameter :: Here = 'elastic (neutronCEstd_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') + 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 @@ -332,10 +353,10 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) class(uncorrelatedReactionCE), pointer :: reac character(100),parameter :: Here =' inelastic (neutronCEstd_class.f90)' - ! Invert inelastic scattering and Get reaction - collDat % MT = self % nuc % invertInelastic(p % 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") + ! Invert inelastic scattering and get reaction + 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") ! Scatter particle if (reac % inCMFrame()) then @@ -402,15 +423,15 @@ subroutine scatterFromFixed(self, p, collDat, reac) integer(shortInt) :: MT ! Read data - MT = collDat % MT + MT = collDat % MT - ! Sample mu , phi and outgoing energy + ! Sample mu, phi and outgoing energy call reac % sampleOut(mu, phi, E_outCM, p % E, p % pRNG) ! 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 call asymptoticInelasticScatter(E_out, mu, E_outCM, collDat % A) diff --git a/CollisionOperator/scatteringKernels_func.f90 b/CollisionOperator/scatteringKernels_func.f90 index 400424fbe..c80df4afb 100644 --- a/CollisionOperator/scatteringKernels_func.f90 +++ b/CollisionOperator/scatteringKernels_func.f90 @@ -12,9 +12,11 @@ module scatteringKernels_func private public :: asymptoticScatter + public :: asymptoticInelasticScatter public :: targetVelocity_constXS public :: targetVelocity_DBRCXS - public :: asymptoticInelasticScatter + public :: relativeEnergy_constXS + public :: dopplerCorrectionFactor private :: sample_x2expx2 private :: sample_x3expx2 @@ -113,7 +115,7 @@ function targetVelocity_constXS(E,dir,A,kT,rand) result (V_t) Y = sqrt(A*E/kT) ! Calculate threshold factor alpha - alpha = 2.0/(Y*sqrt(PI)+2.0) + alpha = TWO / (Y * SQRT_PI + TWO) rejectionLoop: do @@ -125,9 +127,9 @@ function targetVelocity_constXS(E,dir,A,kT,rand) result (V_t) end do rejectionLoop - ! Calculate azimithal angle for target and obtain target direction + ! Calculate azimuthal angle for target and obtain target direction r1 = rand % get() - phi = 2.0 *PI * r1 + phi = TWO * PI * r1 V_t = rotateVector(dir, mu, phi) @@ -164,7 +166,7 @@ function targetVelocity_DBRCXS(nuc, E, dir, A, kT, rand, tempMaj) result (V_t) ! Calculate threshold factor alpha ! In MCNP, alpha is p1 - alpha = 2.0 / (Y * sqrt(PI) + 2.0) + alpha = TWO / (Y * SQRT_PI + TWO) rejectionLoop: do @@ -191,7 +193,7 @@ function targetVelocity_DBRCXS(nuc, E, dir, A, kT, rand, tempMaj) result (V_t) end do rejectionLoop - ! Calculate azimithal angle for target and obtain target direction + ! Calculate azimuthal angle for target and obtain target direction r2 = rand % get() phi = 2.0 * PI * r2 @@ -202,6 +204,73 @@ function targetVelocity_DBRCXS(nuc, E, dir, A, kT, rand, tempMaj) result (V_t) end function targetVelocity_DBRCXS + !! + !! Function that returns a sample of target frame neutron energy using constant XS approximation + !! V_t is a vector. The velocity is scaled by a factor sqrt(Mn/2) where Mn is mass of a neutron + !! so that V_t*V_t=E_t with E_t beeing kinetic energy of a NEUTRON traveling with TARGET VELOCITY + !! (note that it is not a kinetic energy of the target). + !! + function relativeEnergy_constXS(E, A, kT, rand) result(relE) + real(defReal), intent(in) :: E + real(defReal), intent(in) :: A + real(defReal), intent(in) :: kT + class(RNG), intent(inout) :: rand + logical(defBool) :: accept + real(defReal) :: alpha, mu + real(defReal) :: X, Y + real(defReal) :: relV, relE + + ! Calculate neutron Y = beta *V_n + ! beta = sqrt(A*Mn/2kT). Note velocity scaling by sqrt(Mn/2). + Y = sqrt(A*E/kT) + + ! Calculate threshold factor alpha + alpha = TWO / (Y * SQRT_PI + TWO) + + rejectionLoop: do + + ! Sample velocity and calculate angle and acceptance probability + call sample_targetVelocity(X, accept, relV, mu, rand, Y, alpha) + + ! Accept or reject mu + if (accept) exit rejectionLoop + + end do rejectionLoop + + ! Relative energy + relE = (relV**2 * kT / A) + + end function relativeEnergy_constXS + + !! + !! Returns the Doppler Broadening low energy correction factor. When performing Doppler + !! broadening (e.g., TMS), this is multiplied to the cross sections at the base temperature + !! to give the effective cross section. + !! + !! Common notation for this constant is g_E(E, A, kT) or g(v, A, kT). + !! + !! The energy Limits are taken from Serpent 2.1.31 + !! + function dopplerCorrectionFactor(E, A, kT) result(g) + real(defReal), intent(in) :: E + real(defReal), intent(in) :: A + real(defReal), intent(in) :: kT + real(defReal) :: alpha, invAlph + real(defReal) :: g + + alpha = sqrt(A * E / kT) + + invAlph = ONE / alpha + + if (alpha > 250.0_defReal) then + g = ONE + else if (alpha > 2.568_defReal) then + g = ONE + HALF * invAlph * invAlph + else + g = (ONE + HALF * invAlph * invAlph) * erf(alpha) + exp(-alpha * alpha) * invAlph / SQRT_PI + end if + + end function dopplerCorrectionFactor !! !! Helper function to sample x^2 * exp( - x^2) probability distribution @@ -295,10 +364,10 @@ subroutine sample_targetVelocity(X, accept, rel_v, mu, rand, Y, alpha) end if ! Sample polar angle of target velocity wrt. neutron direction - mu = 2.0 * r2 - 1.0; + mu = TWO * r2 - ONE; ! Calculate relative velocity between neutron and target - rel_v = sqrt(Y * Y + X * X - 2.0 * X * Y * mu) + rel_v = sqrt(Y * Y + X * X - TWO * X * Y * mu) ! Calculate Acceptance Propability P_acc = rel_v / (Y + X) @@ -308,5 +377,4 @@ subroutine sample_targetVelocity(X, accept, rel_v, mu, rand, Y, alpha) end subroutine sample_targetVelocity - end module scatteringKernels_func diff --git a/InputFiles/Benchmarks/BEAVRS/BEAVRS b/InputFiles/Benchmarks/BEAVRS/BEAVRS index 3c7cb89de..54ac5f845 100644 --- a/InputFiles/Benchmarks/BEAVRS/BEAVRS +++ b/InputFiles/Benchmarks/BEAVRS/BEAVRS @@ -1954,7 +1954,7 @@ nuclearData { // but are not available in the JEFF-3.11 data library Air { - temp 566; + !temp 566; composition { 18036.06 7.8730E-09; 18038.06 1.4844E-09; @@ -1972,27 +1972,27 @@ nuclearData { SS304 { temp 566; composition { - 24050.06 7.6778E-04; - 24052.06 1.4806E-02; - 24053.06 1.6789E-03; - 24054.06 4.1791E-04; - 26054.06 3.4620E-03; - 26056.06 5.4345E-02; - 26057.06 1.2551E-03; - 26058.06 1.6703E-04; - 25055.06 1.7604E-03; - 28058.06 5.6089E-03; - 28060.06 2.1605E-03; - 28061.06 9.3917E-05; - 28062.06 2.9945E-04; - 28064.06 7.6261E-05; - 14028.06 9.5281E-04; - 14029.06 4.8381E-05; - 14030.06 3.1893E-05; } + 24050.03 7.6778E-04; + 24052.03 1.4806E-02; + 24053.03 1.6789E-03; + 24054.03 4.1791E-04; + 26054.03 3.4620E-03; + 26056.03 5.4345E-02; + 26057.03 1.2551E-03; + 26058.03 1.6703E-04; + 25055.03 1.7604E-03; + 28058.03 5.6089E-03; + 28060.03 2.1605E-03; + 28061.03 9.3917E-05; + 28062.03 2.9945E-04; + 28064.03 7.6261E-05; + 14028.03 9.5281E-04; + 14029.03 4.8381E-05; + 14030.03 3.1893E-05; } } Helium { - temp 566; + !temp 566; composition { 2003.06 4.8089E-10; 2004.06 2.4044E-04; } @@ -2001,20 +2001,20 @@ nuclearData { BorosilicateGlass { temp 566; composition { - 13027.06 1.7352E-03; - 5010.06 9.6506E-04; - 5011.06 3.9189E-03; - 8016.06 4.6514E-02; - 8017.06 1.7671E-05; - //8018.06 9.3268E-05; - 14028.06 1.6926E-02; - 14029.06 8.5944E-04; - 14030.06 5.6654E-04; } + 13027.03 1.7352E-03; + 5010.03 9.6506E-04; + 5011.03 3.9189E-03; + 8016.03 4.6514E-02; + 8017.03 1.7671E-05; + //8018.03 9.3268E-05; + 14028.03 1.6926E-02; + 14029.03 8.5944E-04; + 14030.03 5.6654E-04; } } Water { - temp 566; - moder {1001.06 lwj3.11; } + !temp 566; + moder {1001.03 lwj3.11; } composition { 5010.06 7.9714E-06; 5011.06 3.2247E-05; @@ -2029,58 +2029,58 @@ nuclearData { Zircaloy { temp 566; composition { - 24050.06 3.2962E-06; - 24052.06 6.3564E-05; - 24053.06 7.2076E-06; - 24054.06 1.7941E-06; - 26054.06 8.6698E-06; - 26056.06 1.3610E-04; - 26057.06 3.1431E-06; - 26058.06 4.1829E-07; - 8016.06 3.0744E-04; - 8017.06 1.1680E-07; - //8018.06 6.1648E-07; - 50112.06 4.6735E-06; - 50114.06 3.1799E-06; - 50115.06 1.6381E-06; - 50116.06 7.0055E-05; - 50117.06 3.7003E-05; - 50118.06 1.1669E-04; - 50119.06 4.1387E-05; - 50120.06 1.5697E-04; - 50122.06 2.2308E-05; - 50124.06 2.7897E-05; - 40090.06 2.1828E-02; - 40091.06 4.7601E-03; - 40092.06 7.2759E-03; - 40094.06 7.3734E-03; - 40096.06 1.1879E-03; } + 24050.03 3.2962E-06; + 24052.03 6.3564E-05; + 24053.03 7.2076E-06; + 24054.03 1.7941E-06; + 26054.03 8.6698E-06; + 26056.03 1.3610E-04; + 26057.03 3.1431E-06; + 26058.03 4.1829E-07; + 8016.03 3.0744E-04; + 8017.03 1.1680E-07; + //8018.03 6.1648E-07; + 50112.03 4.6735E-06; + 50114.03 3.1799E-06; + 50115.03 1.6381E-06; + 50116.03 7.0055E-05; + 50117.03 3.7003E-05; + 50118.03 1.1669E-04; + 50119.03 4.1387E-05; + 50120.03 1.5697E-04; + 50122.03 2.2308E-05; + 50124.03 2.7897E-05; + 40090.03 2.1828E-02; + 40091.03 4.7601E-03; + 40092.03 7.2759E-03; + 40094.03 7.3734E-03; + 40096.03 1.1879E-03; } } Inconel{ temp 566; composition { - 24050.06 7.8239E-04; - 24052.06 1.5088E-02; - 24053.06 1.7108E-03; - 24054.06 4.2586E-04; - 26054.06 1.4797E-03; - 26056.06 2.3229E-02; - 26057.06 5.3645E-04; - 26058.06 7.1392E-05; - 25055.06 7.8201E-04; - 28058.06 2.9320E-02; - 28060.06 1.1294E-02; - 28061.06 4.9094E-04; - 28062.06 1.5653E-03; - 28064.06 3.9864E-04; - 14028.06 5.6757E-04; - 14029.06 2.8820E-05; - 14030.06 1.8998E-05; } + 24050.03 7.8239E-04; + 24052.03 1.5088E-02; + 24053.03 1.7108E-03; + 24054.03 4.2586E-04; + 26054.03 1.4797E-03; + 26056.03 2.3229E-02; + 26057.03 5.3645E-04; + 26058.03 7.1392E-05; + 25055.03 7.8201E-04; + 28058.03 2.9320E-02; + 28060.03 1.1294E-02; + 28061.03 4.9094E-04; + 28062.03 1.5653E-03; + 28064.03 3.9864E-04; + 14028.03 5.6757E-04; + 14029.03 2.8820E-05; + 14030.03 1.8998E-05; } } B4C{ - temp 566; + !temp 566; composition { 5010.06 1.5206E-02; 5011.06 6.1514E-02; @@ -2092,75 +2092,75 @@ nuclearData { Ag-In-Cd{ temp 566; composition { - 47107.06 2.3523E-02; - 47109.06 2.1854E-02; - 48106.06 3.3882E-05; - 48108.06 2.4166E-05; - 48110.06 3.3936E-04; - 48111.06 3.4821E-04; - 48112.06 6.5611E-04; - 48113.06 3.3275E-04; - 48114.06 7.8252E-04; - 48116.06 2.0443E-04; - 49113.06 3.4219E-04; - 49115.06 7.6511E-03; } + 47107.03 2.3523E-02; + 47109.03 2.1854E-02; + 48106.03 3.3882E-05; + 48108.03 2.4166E-05; + 48110.03 3.3936E-04; + 48111.03 3.4821E-04; + 48112.03 6.5611E-04; + 48113.03 3.3275E-04; + 48114.03 7.8252E-04; + 48116.03 2.0443E-04; + 49113.03 3.4219E-04; + 49115.03 7.6511E-03; } } //mox fuel used UO2-16 { temp 566; composition { - 8016.06 4.5897E-02; - 8017.06 1.7436E-05; - //8018.06 9.2032E-05; - 92234.06 3.0131E-06; - 92235.06 3.7503E-04; - 92238.06 2.2625E-02;} + 8016.03 4.5897E-02; + 8017.03 1.7436E-05; + //8018.03 9.2032E-05; + 92234.03 3.0131E-06; + 92235.03 3.7503E-04; + 92238.03 2.2625E-02;} } UO2-24 { temp 566; composition { - 8016.06 4.5830E-02; - 8017.06 1.7411E-05; - //8018.06 9.1898E-05; - 92234.06 4.4842E-06; - 92235.06 5.5814E-04; - 92238.06 2.2407E-02;} + 8016.03 4.5830E-02; + 8017.03 1.7411E-05; + //8018.03 9.1898E-05; + 92234.03 4.4842E-06; + 92235.03 5.5814E-04; + 92238.03 2.2407E-02;} } UO2-31 { temp 566; composition { - 8016.06 4.5853E-02; - 8017.06 1.7420E-05; - //8018.06 9.1942E-05; - 92234.06 5.7987E-06; - 92235.06 7.2175E-04; - 92238.06 2.2253E-02;} + 8016.03 4.5853E-02; + 8017.03 1.7420E-05; + //8018.03 9.1942E-05; + 92234.03 5.7987E-06; + 92235.03 7.2175E-04; + 92238.03 2.2253E-02;} } UO2-32 { temp 566; composition { - 8016.06 4.6029E-02; - 8017.06 1.7487E-05; - //8018.06 9.2296E-05; - 92234.06 5.9959E-06; - 92235.06 7.4630E-04; - 92238.06 2.2317E-02; + 8016.03 4.6029E-02; + 8017.03 1.7487E-05; + //8018.03 9.2296E-05; + 92234.03 5.9959E-06; + 92235.03 7.4630E-04; + 92238.03 2.2317E-02; } } UO2-34 { temp 566; composition { - 8016.06 4.6110E-02; - 8017.06 1.7517E-05; - //8018.06 9.2459E-05; - 92234.06 6.4018E-06; - 92235.06 7.9681E-04; - 92238.06 2.2307E-02;} + 8016.03 4.6110E-02; + 8017.03 1.7517E-05; + //8018.03 9.2459E-05; + 92234.03 6.4018E-06; + 92235.03 7.9681E-04; + 92238.03 2.2307E-02;} } // vanadium51 was stated twice in carbonsteel below @@ -2168,92 +2168,92 @@ nuclearData { CarbonSteel { temp 566; composition { - 13027.06 4.3523E-05; - 5010.06 2.5833E-06; - 5011.06 1.0450E-05; - 6012.06 1.0442E-03; - //6013.06 1.1697E-05 ; - 20040.06 1.7043E-05; - 20042.06 1.1375E-07; - 20043.06 2.3734E-08; - 20044.06 3.6673E-07; - 20046.06 7.0322E-10; - 20048.06 3.2875E-08; - 24050.06 1.3738E-05; - 24052.06 2.6493E-04; - 24053.06 3.0041E-05; - 24054.06 7.4778E-06; - 29063.06 1.0223E-04; - 29065.06 4.5608E-05; - 26054.06 4.7437E-03; - 26056.06 7.4465E-02; - 26057.06 1.7197E-03; - 26058.06 2.2886E-04; - 25055.06 6.4126E-04; - 42100.06 2.9814E-05; - 42092.06 4.4822E-05; - 42094.06 2.8110E-05; - 42095.06 4.8567E-05; - 42096.06 5.1015E-05; - 42097.06 2.9319E-05; - 42098.06 7.4327E-05; - 41093.06 5.0559E-06; - 28058.06 4.0862E-04; - 28060.06 1.5740E-04; - 28061.06 6.8420E-06; - 28062.06 2.1815E-05; - 28064.06 5.5557E-06; - 15031.06 3.7913E-05; - 16032.06 3.4808E-05; - 16033.06 2.7420E-07; - 16034.06 1.5368E-06; - 16036.06 5.3398E-09; - 14028.06 6.1702E-04; - 14029.06 3.1330E-05; - 14030.06 2.0653E-05; - 22046.06 1.2144E-06; - 22047.06 1.0952E-06; - 22048.06 1.0851E-05; - 22049.06 7.9634E-07; - 22050.06 7.6249E-07; - //23050.06 1.1526E-07; - 23051.06 4.5989E-05; + 13027.03 4.3523E-05; + 5010.03 2.5833E-06; + 5011.03 1.0450E-05; + 6012.03 1.0442E-03; + //6013.03 1.1697E-05 ; + 20040.03 1.7043E-05; + 20042.03 1.1375E-07; + 20043.03 2.3734E-08; + 20044.03 3.6673E-07; + 20046.03 7.0322E-10; + 20048.03 3.2875E-08; + 24050.03 1.3738E-05; + 24052.03 2.6493E-04; + 24053.03 3.0041E-05; + 24054.03 7.4778E-06; + 29063.03 1.0223E-04; + 29065.03 4.5608E-05; + 26054.03 4.7437E-03; + 26056.03 7.4465E-02; + 26057.03 1.7197E-03; + 26058.03 2.2886E-04; + 25055.03 6.4126E-04; + 42100.03 2.9814E-05; + 42092.03 4.4822E-05; + 42094.03 2.8110E-05; + 42095.03 4.8567E-05; + 42096.03 5.1015E-05; + 42097.03 2.9319E-05; + 42098.03 7.4327E-05; + 41093.03 5.0559E-06; + 28058.03 4.0862E-04; + 28060.03 1.5740E-04; + 28061.03 6.8420E-06; + 28062.03 2.1815E-05; + 28064.03 5.5557E-06; + 15031.03 3.7913E-05; + 16032.03 3.4808E-05; + 16033.03 2.7420E-07; + 16034.03 1.5368E-06; + 16036.03 5.3398E-09; + 14028.03 6.1702E-04; + 14029.03 3.1330E-05; + 14030.03 2.0653E-05; + 22046.03 1.2144E-06; + 22047.03 1.0952E-06; + 22048.03 1.0851E-05; + 22049.03 7.9634E-07; + 22050.03 7.6249E-07; + //23050.03 1.1526E-07; + 23051.03 4.5989E-05; } } SupportPlateSS { temp 566; composition { - 24050.06 3.5223E-04; - 24052.06 6.7924E-03; - 24053.06 7.7020E-04; - 24054.06 1.9172E-04; - 26054.06 1.5882E-03; - 26056.06 2.4931E-02; - 26057.06 5.7578E-04; - 26058.06 7.6625E-05; - 25055.06 8.0762E-04; - 28058.06 2.5731E-03; - 28060.06 9.9117E-04; - 28061.06 4.3085E-05; - 28062.06 1.3738E-04; - 28064.06 3.4985E-05; - 14028.06 4.3711E-04; - 14029.06 2.2195E-05; - 14030.06 1.4631E-05;} + 24050.03 3.5223E-04; + 24052.03 6.7924E-03; + 24053.03 7.7020E-04; + 24054.03 1.9172E-04; + 26054.03 1.5882E-03; + 26056.03 2.4931E-02; + 26057.03 5.7578E-04; + 26058.03 7.6625E-05; + 25055.03 8.0762E-04; + 28058.03 2.5731E-03; + 28060.03 9.9117E-04; + 28061.03 4.3085E-05; + 28062.03 1.3738E-04; + 28064.03 3.4985E-05; + 14028.03 4.3711E-04; + 14029.03 2.2195E-05; + 14030.03 1.4631E-05;} } SupportPlateBW { temp 566; - moder {1001.06 lwj3.11; } + moder {1001.03 lwj3.11; } composition { - 5010.06 1.0559E-05; - 5011.06 4.2716E-05; - 1001.06 6.5512E-02; - 1002.06 1.0204E-05; - 8016.06 3.2683E-02; - 8017.06 1.2416E-05; - //8018.06 6.5535E-05; + 5010.03 1.0559E-05; + 5011.03 4.2716E-05; + 1001.03 6.5512E-02; + 1002.03 1.0204E-05; + 8016.03 3.2683E-02; + 8017.03 1.2416E-05; + //8018.03 6.5535E-05; } } diff --git a/InputFiles/Benchmarks/BEAVRS/BEAVRS2D b/InputFiles/Benchmarks/BEAVRS/BEAVRS2D index 5a9bac1d1..441d43e6d 100644 --- a/InputFiles/Benchmarks/BEAVRS/BEAVRS2D +++ b/InputFiles/Benchmarks/BEAVRS/BEAVRS2D @@ -1273,7 +1273,7 @@ nuclearData { ! These are commented out. Air { - temp 566; + !temp 566; composition { 18036.06 7.8730E-09; 18038.06 1.4844E-09; @@ -1291,27 +1291,27 @@ nuclearData { SS304 { temp 566; composition { - 24050.06 7.6778E-04; - 24052.06 1.4806E-02; - 24053.06 1.6789E-03; - 24054.06 4.1791E-04; - 26054.06 3.4620E-03; - 26056.06 5.4345E-02; - 26057.06 1.2551E-03; - 26058.06 1.6703E-04; - 25055.06 1.7604E-03; - 28058.06 5.6089E-03; - 28060.06 2.1605E-03; - 28061.06 9.3917E-05; - 28062.06 2.9945E-04; - 28064.06 7.6261E-05; - 14028.06 9.5281E-04; - 14029.06 4.8381E-05; - 14030.06 3.1893E-05; } + 24050.03 7.6778E-04; + 24052.03 1.4806E-02; + 24053.03 1.6789E-03; + 24054.03 4.1791E-04; + 26054.03 3.4620E-03; + 26056.03 5.4345E-02; + 26057.03 1.2551E-03; + 26058.03 1.6703E-04; + 25055.03 1.7604E-03; + 28058.03 5.6089E-03; + 28060.03 2.1605E-03; + 28061.03 9.3917E-05; + 28062.03 2.9945E-04; + 28064.03 7.6261E-05; + 14028.03 9.5281E-04; + 14029.03 4.8381E-05; + 14030.03 3.1893E-05; } } Helium { - temp 566; + !temp 566; composition { 2003.06 4.8089E-10; 2004.06 2.4044E-04; } @@ -1320,20 +1320,20 @@ nuclearData { BorosilicateGlass { temp 566; composition { - 13027.06 1.7352E-03; - 5010.06 9.6506E-04; - 5011.06 3.9189E-03; - 8016.06 4.6514E-02; - 8017.06 1.7671E-05; - //8018.06 9.3268E-05; - 14028.06 1.6926E-02; - 14029.06 8.5944E-04; - 14030.06 5.6654E-04; } + 13027.03 1.7352E-03; + 5010.03 9.6506E-04; + 5011.03 3.9189E-03; + 8016.03 4.6514E-02; + 8017.03 1.7671E-05; + //8018.03 9.3268E-05; + 14028.03 1.6926E-02; + 14029.03 8.5944E-04; + 14030.03 5.6654E-04; } } Water { - temp 566; - moder {1001.06 lwj3.11; } + !temp 566; + moder {1001.03 lwj3.11; } composition { 5010.06 7.9714E-06; 5011.06 3.2247E-05; @@ -1348,58 +1348,58 @@ nuclearData { Zircaloy { temp 566; composition { - 24050.06 3.2962E-06; - 24052.06 6.3564E-05; - 24053.06 7.2076E-06; - 24054.06 1.7941E-06; - 26054.06 8.6698E-06; - 26056.06 1.3610E-04; - 26057.06 3.1431E-06; - 26058.06 4.1829E-07; - 8016.06 3.0744E-04; - 8017.06 1.1680E-07; - //8018.06 6.1648E-07; - 50112.06 4.6735E-06; - 50114.06 3.1799E-06; - 50115.06 1.6381E-06; - 50116.06 7.0055E-05; - 50117.06 3.7003E-05; - 50118.06 1.1669E-04; - 50119.06 4.1387E-05; - 50120.06 1.5697E-04; - 50122.06 2.2308E-05; - 50124.06 2.7897E-05; - 40090.06 2.1828E-02; - 40091.06 4.7601E-03; - 40092.06 7.2759E-03; - 40094.06 7.3734E-03; - 40096.06 1.1879E-03; } + 24050.03 3.2962E-06; + 24052.03 6.3564E-05; + 24053.03 7.2076E-06; + 24054.03 1.7941E-06; + 26054.03 8.6698E-06; + 26056.03 1.3610E-04; + 26057.03 3.1431E-06; + 26058.03 4.1829E-07; + 8016.03 3.0744E-04; + 8017.03 1.1680E-07; + //8018.03 6.1648E-07; + 50112.03 4.6735E-06; + 50114.03 3.1799E-06; + 50115.03 1.6381E-06; + 50116.03 7.0055E-05; + 50117.03 3.7003E-05; + 50118.03 1.1669E-04; + 50119.03 4.1387E-05; + 50120.03 1.5697E-04; + 50122.03 2.2308E-05; + 50124.03 2.7897E-05; + 40090.03 2.1828E-02; + 40091.03 4.7601E-03; + 40092.03 7.2759E-03; + 40094.03 7.3734E-03; + 40096.03 1.1879E-03; } } Inconel{ temp 566; composition { - 24050.06 7.8239E-04; - 24052.06 1.5088E-02; - 24053.06 1.7108E-03; - 24054.06 4.2586E-04; - 26054.06 1.4797E-03; - 26056.06 2.3229E-02; - 26057.06 5.3645E-04; - 26058.06 7.1392E-05; - 25055.06 7.8201E-04; - 28058.06 2.9320E-02; - 28060.06 1.1294E-02; - 28061.06 4.9094E-04; - 28062.06 1.5653E-03; - 28064.06 3.9864E-04; - 14028.06 5.6757E-04; - 14029.06 2.8820E-05; - 14030.06 1.8998E-05; } + 24050.03 7.8239E-04; + 24052.03 1.5088E-02; + 24053.03 1.7108E-03; + 24054.03 4.2586E-04; + 26054.03 1.4797E-03; + 26056.03 2.3229E-02; + 26057.03 5.3645E-04; + 26058.03 7.1392E-05; + 25055.03 7.8201E-04; + 28058.03 2.9320E-02; + 28060.03 1.1294E-02; + 28061.03 4.9094E-04; + 28062.03 1.5653E-03; + 28064.03 3.9864E-04; + 14028.03 5.6757E-04; + 14029.03 2.8820E-05; + 14030.03 1.8998E-05; } } B4C{ - temp 566; + !temp 566; composition { 5010.06 1.5206E-02; 5011.06 6.1514E-02; @@ -1411,75 +1411,75 @@ nuclearData { Ag-In-Cd{ temp 566; composition { - 47107.06 2.3523E-02; - 47109.06 2.1854E-02; - 48106.06 3.3882E-05; - 48108.06 2.4166E-05; - 48110.06 3.3936E-04; - 48111.06 3.4821E-04; - 48112.06 6.5611E-04; - 48113.06 3.3275E-04; - 48114.06 7.8252E-04; - 48116.06 2.0443E-04; - 49113.06 3.4219E-04; - 49115.06 7.6511E-03; } + 47107.03 2.3523E-02; + 47109.03 2.1854E-02; + 48106.03 3.3882E-05; + 48108.03 2.4166E-05; + 48110.03 3.3936E-04; + 48111.03 3.4821E-04; + 48112.03 6.5611E-04; + 48113.03 3.3275E-04; + 48114.03 7.8252E-04; + 48116.03 2.0443E-04; + 49113.03 3.4219E-04; + 49115.03 7.6511E-03; } } //mox fuel used UO2-16 { temp 566; composition { - 8016.06 4.5897E-02; - 8017.06 1.7436E-05; - //8018.06 9.2032E-05; - 92234.06 3.0131E-06; - 92235.06 3.7503E-04; - 92238.06 2.2625E-02;} + 8016.03 4.5897E-02; + 8017.03 1.7436E-05; + //8018.03 9.2032E-05; + 92234.03 3.0131E-06; + 92235.03 3.7503E-04; + 92238.03 2.2625E-02;} } UO2-24 { temp 566; composition { - 8016.06 4.5830E-02; - 8017.06 1.7411E-05; - //8018.06 9.1898E-05; - 92234.06 4.4842E-06; - 92235.06 5.5814E-04; - 92238.06 2.2407E-02;} + 8016.03 4.5830E-02; + 8017.03 1.7411E-05; + //8018.03 9.1898E-05; + 92234.03 4.4842E-06; + 92235.03 5.5814E-04; + 92238.03 2.2407E-02;} } UO2-31 { temp 566; composition { - 8016.06 4.5853E-02; - 8017.06 1.7420E-05; - //8018.06 9.1942E-05; - 92234.06 5.7987E-06; - 92235.06 7.2175E-04; - 92238.06 2.2253E-02;} + 8016.03 4.5853E-02; + 8017.03 1.7420E-05; + //8018.03 9.1942E-05; + 92234.03 5.7987E-06; + 92235.03 7.2175E-04; + 92238.03 2.2253E-02;} } UO2-32 { temp 566; composition { - 8016.06 4.6029E-02; - 8017.06 1.7487E-05; - //8018.06 9.2296E-05; - 92234.06 5.9959E-06; - 92235.06 7.4630E-04; - 92238.06 2.2317E-02; + 8016.03 4.6029E-02; + 8017.03 1.7487E-05; + //8018.03 9.2296E-05; + 92234.03 5.9959E-06; + 92235.03 7.4630E-04; + 92238.03 2.2317E-02; } } UO2-34 { temp 566; composition { - 8016.06 4.6110E-02; - 8017.06 1.7517E-05; - //8018.06 9.2459E-05; - 92234.06 6.4018E-06; - 92235.06 7.9681E-04; - 92238.06 2.2307E-02;} + 8016.03 4.6110E-02; + 8017.03 1.7517E-05; + //8018.03 9.2459E-05; + 92234.03 6.4018E-06; + 92235.03 7.9681E-04; + 92238.03 2.2307E-02;} } // vanadium51 was stated twice in carbonsteel below @@ -1487,92 +1487,92 @@ nuclearData { CarbonSteel { temp 566; composition { - 13027.06 4.3523E-05; - 5010.06 2.5833E-06; - 5011.06 1.0450E-05; - 6012.06 1.0442E-03; - //6013.06 1.1697E-05 ; - 20040.06 1.7043E-05; - 20042.06 1.1375E-07; - 20043.06 2.3734E-08; - 20044.06 3.6673E-07; - 20046.06 7.0322E-10; - 20048.06 3.2875E-08; - 24050.06 1.3738E-05; - 24052.06 2.6493E-04; - 24053.06 3.0041E-05; - 24054.06 7.4778E-06; - 29063.06 1.0223E-04; - 29065.06 4.5608E-05; - 26054.06 4.7437E-03; - 26056.06 7.4465E-02; - 26057.06 1.7197E-03; - 26058.06 2.2886E-04; - 25055.06 6.4126E-04; - 42100.06 2.9814E-05; - 42092.06 4.4822E-05; - 42094.06 2.8110E-05; - 42095.06 4.8567E-05; - 42096.06 5.1015E-05; - 42097.06 2.9319E-05; - 42098.06 7.4327E-05; - 41093.06 5.0559E-06; - 28058.06 4.0862E-04; - 28060.06 1.5740E-04; - 28061.06 6.8420E-06; - 28062.06 2.1815E-05; - 28064.06 5.5557E-06; - 15031.06 3.7913E-05; - 16032.06 3.4808E-05; - 16033.06 2.7420E-07; - 16034.06 1.5368E-06; - 16036.06 5.3398E-09; - 14028.06 6.1702E-04; - 14029.06 3.1330E-05; - 14030.06 2.0653E-05; - 22046.06 1.2144E-06; - 22047.06 1.0952E-06; - 22048.06 1.0851E-05; - 22049.06 7.9634E-07; - 22050.06 7.6249E-07; - //23050.06 1.1526E-07; - 23051.06 4.5989E-05; + 13027.03 4.3523E-05; + 5010.03 2.5833E-06; + 5011.03 1.0450E-05; + 6012.03 1.0442E-03; + //6013.03 1.1697E-05 ; + 20040.03 1.7043E-05; + 20042.03 1.1375E-07; + 20043.03 2.3734E-08; + 20044.03 3.6673E-07; + 20046.03 7.0322E-10; + 20048.03 3.2875E-08; + 24050.03 1.3738E-05; + 24052.03 2.6493E-04; + 24053.03 3.0041E-05; + 24054.03 7.4778E-06; + 29063.03 1.0223E-04; + 29065.03 4.5608E-05; + 26054.03 4.7437E-03; + 26056.03 7.4465E-02; + 26057.03 1.7197E-03; + 26058.03 2.2886E-04; + 25055.03 6.4126E-04; + 42100.03 2.9814E-05; + 42092.03 4.4822E-05; + 42094.03 2.8110E-05; + 42095.03 4.8567E-05; + 42096.03 5.1015E-05; + 42097.03 2.9319E-05; + 42098.03 7.4327E-05; + 41093.03 5.0559E-06; + 28058.03 4.0862E-04; + 28060.03 1.5740E-04; + 28061.03 6.8420E-06; + 28062.03 2.1815E-05; + 28064.03 5.5557E-06; + 15031.03 3.7913E-05; + 16032.03 3.4808E-05; + 16033.03 2.7420E-07; + 16034.03 1.5368E-06; + 16036.03 5.3398E-09; + 14028.03 6.1702E-04; + 14029.03 3.1330E-05; + 14030.03 2.0653E-05; + 22046.03 1.2144E-06; + 22047.03 1.0952E-06; + 22048.03 1.0851E-05; + 22049.03 7.9634E-07; + 22050.03 7.6249E-07; + //23050.03 1.1526E-07; + 23051.03 4.5989E-05; } } SupportPlateSS { temp 566; composition { - 24050.06 3.5223E-04; - 24052.06 6.7924E-03; - 24053.06 7.7020E-04; - 24054.06 1.9172E-04; - 26054.06 1.5882E-03; - 26056.06 2.4931E-02; - 26057.06 5.7578E-04; - 26058.06 7.6625E-05; - 25055.06 8.0762E-04; - 28058.06 2.5731E-03; - 28060.06 9.9117E-04; - 28061.06 4.3085E-05; - 28062.06 1.3738E-04; - 28064.06 3.4985E-05; - 14028.06 4.3711E-04; - 14029.06 2.2195E-05; - 14030.06 1.4631E-05;} + 24050.03 3.5223E-04; + 24052.03 6.7924E-03; + 24053.03 7.7020E-04; + 24054.03 1.9172E-04; + 26054.03 1.5882E-03; + 26056.03 2.4931E-02; + 26057.03 5.7578E-04; + 26058.03 7.6625E-05; + 25055.03 8.0762E-04; + 28058.03 2.5731E-03; + 28060.03 9.9117E-04; + 28061.03 4.3085E-05; + 28062.03 1.3738E-04; + 28064.03 3.4985E-05; + 14028.03 4.3711E-04; + 14029.03 2.2195E-05; + 14030.03 1.4631E-05;} } SupportPlateBW { temp 566; - moder {1001.06 lwj3.11; } + moder {1001.03 lwj3.11; } composition { - 5010.06 1.0559E-05; - 5011.06 4.2716E-05; - 1001.06 6.5512E-02; - 1002.06 1.0204E-05; - 8016.06 3.2683E-02; - 8017.06 1.2416E-05; - //8018.06 6.5535E-05; + 5010.03 1.0559E-05; + 5011.03 4.2716E-05; + 1001.03 6.5512E-02; + 1002.03 1.0204E-05; + 8016.03 3.2683E-02; + 8017.03 1.2416E-05; + //8018.03 6.5535E-05; } } diff --git a/InputFiles/Benchmarks/HoogenboomMartin/HMcore b/InputFiles/Benchmarks/HoogenboomMartin/HMcore index b57d995ca..1ab1e93fd 100644 --- a/InputFiles/Benchmarks/HoogenboomMartin/HMcore +++ b/InputFiles/Benchmarks/HoogenboomMartin/HMcore @@ -289,7 +289,6 @@ materials { UOXfuel { - temp 16486; composition { 92234.03 4.9476e-6; 92235.03 4.8218e-4; @@ -328,7 +327,6 @@ UOXfuel { } Zrcladding { - temp 31628; composition { 40090.03 0.019578; 40091.03 0.0042689; @@ -338,7 +336,6 @@ Zrcladding { } coldwater { - temp 45286; moder { 1001.03 lwj3.00; } composition { 1001.03 0.04938; @@ -350,7 +347,6 @@ coldwater { hotwater { - temp 46865; moder { 1001.03 lwj3.00; } composition { 1001.03 0.04404; @@ -363,7 +359,6 @@ hotwater { reactorvessel { - temp 51694; composition { 26054.03 4.79006E-03; 26056.03 7.5184E-02; @@ -397,7 +392,6 @@ reactorvessel { reflectorbottom { - temp 56982; moder { 1001.03 lwj3.00; } composition { 1001.03 2.4886E-02; @@ -426,7 +420,6 @@ reflectorbottom { reflectortop { - temp 54357; moder { 1001.03 lwj3.00; } composition { 1001.03 2.2196E-02; @@ -455,7 +448,6 @@ reflectortop { coreplate { - temp 54687; moder { 1001.03 lwj3.00; } composition { 1001.03 4.9773E-03; @@ -483,7 +475,6 @@ coreplate { } nozzlebottom { - temp 87648; moder { 1001.03 lwj3.00; } composition { 1001.03 3.733E-2; @@ -512,7 +503,6 @@ nozzlebottom { nozzletop { - temp 87064; moder { 1001.03 lwj3.00; } composition { 1001.03 3.7733E-2; @@ -540,7 +530,6 @@ nozzletop { } FAbottom { - temp 45034; moder { 1001.03 lwj3.00; } composition { 1001.03 2.986E-02; @@ -556,7 +545,6 @@ FAbottom { } FAtop { - temp 43053; moder { 1001.03 lwj3.00; } composition { 1001.03 3.107E-02; diff --git a/InputFiles/Benchmarks/MCNP/BIGTEN b/InputFiles/Benchmarks/MCNP/BIGTEN index 777b15685..3f7aaba13 100644 --- a/InputFiles/Benchmarks/MCNP/BIGTEN +++ b/InputFiles/Benchmarks/MCNP/BIGTEN @@ -74,7 +74,7 @@ nuclearData { materials { core { - temp 293; + !temp 293; composition { 92234.03 4.8416E-5; 92235.03 4.8151E-3; @@ -83,7 +83,7 @@ nuclearData { } } refl { - temp 293; + !temp 293; composition { 92234.03 2.8672E-7; 92235.03 1.0058E-4; diff --git a/InputFiles/Benchmarks/MCNP/Falstaff b/InputFiles/Benchmarks/MCNP/Falstaff index bade833ed..b41a6f322 100644 --- a/InputFiles/Benchmarks/MCNP/Falstaff +++ b/InputFiles/Benchmarks/MCNP/Falstaff @@ -74,7 +74,7 @@ nuclearData { materials { solution { - temp 293; + !temp 293; composition { 92232.03 0.000000045608; 92233.03 0.0022379; @@ -88,7 +88,7 @@ materials { moder {01001.03 hh2o.04; } } steel { - temp 293; + !temp 293; composition { 26000.03 0.061248; 24000.03 0.016678; @@ -96,7 +96,7 @@ materials { } } reflector { - temp 293; + !temp 293; composition { 04009.03 0.12161; } diff --git a/InputFiles/Benchmarks/MCNP/IEUMF04 b/InputFiles/Benchmarks/MCNP/IEUMF04 index 554d27119..7f035fca6 100644 --- a/InputFiles/Benchmarks/MCNP/IEUMF04 +++ b/InputFiles/Benchmarks/MCNP/IEUMF04 @@ -88,7 +88,7 @@ nuclearData { layer1 { - temp 293; + !temp 293; composition { 92234.03 1.5926E-4; 92235.03 1.7443E-2; @@ -99,7 +99,7 @@ nuclearData { } } layer2 { - temp 293; + !temp 293; composition { 92234.03 1.5878E-4; 92235.03 1.7415E-2; @@ -110,7 +110,7 @@ nuclearData { } } layer3 { - temp 293; + !temp 293; composition { 92234.03 1.5803E-4; 92235.03 1.7418E-2; @@ -121,7 +121,7 @@ nuclearData { } } layer4 { - temp 293; + !temp 293; composition { 92234.03 1.3423E-4; 92235.03 1.7356E-2; @@ -132,7 +132,7 @@ nuclearData { } } layer5 { - temp 293; + !temp 293; composition { 92234.03 1.6281E-4; 92235.03 1.7417E-2; @@ -143,7 +143,7 @@ nuclearData { } } layer6 { - temp 293; + !temp 293; composition { 92234.03 1.7609E-4; 92235.03 1.7326E-2; @@ -154,7 +154,7 @@ nuclearData { } } layer7 { - temp 293; + !temp 293; composition { 92234.03 1.5156E-4; 92235.03 1.7266E-2; @@ -165,7 +165,7 @@ nuclearData { } } refl { - temp 293; + !temp 293; composition { 06012.03 7.7716E-2; } diff --git a/InputFiles/Benchmarks/MCNP/Jezebel b/InputFiles/Benchmarks/MCNP/Jezebel index 03ec450c0..bedf3ce60 100644 --- a/InputFiles/Benchmarks/MCNP/Jezebel +++ b/InputFiles/Benchmarks/MCNP/Jezebel @@ -76,7 +76,7 @@ materials { fuel { - temp 293; + !temp 293; composition { 94239.03 3.7047E-2; 94240.03 1.7512E-3; diff --git a/InputFiles/Benchmarks/MCNP/Jezebel233 b/InputFiles/Benchmarks/MCNP/Jezebel233 index aa990732a..d6b20c5d8 100644 --- a/InputFiles/Benchmarks/MCNP/Jezebel233 +++ b/InputFiles/Benchmarks/MCNP/Jezebel233 @@ -73,7 +73,7 @@ materials { fuel { - temp 293; + !temp 293; composition { 92233.03 0.046712; 92234.03 0.00059026; diff --git a/InputFiles/Benchmarks/MCNP/Jezebel240 b/InputFiles/Benchmarks/MCNP/Jezebel240 index 600e6761f..d3999f596 100644 --- a/InputFiles/Benchmarks/MCNP/Jezebel240 +++ b/InputFiles/Benchmarks/MCNP/Jezebel240 @@ -76,7 +76,7 @@ materials { fuel { - temp 293; + !temp 293; composition { 94239.03 2.9934E-2; 94240.03 7.8754E-3; diff --git a/InputFiles/Benchmarks/MCNP/LEUST02 b/InputFiles/Benchmarks/MCNP/LEUST02 index 66121dbfd..8ed86030b 100644 --- a/InputFiles/Benchmarks/MCNP/LEUST02 +++ b/InputFiles/Benchmarks/MCNP/LEUST02 @@ -72,7 +72,7 @@ nuclearData { solution { - temp 293; + !temp 293; composition { 92234.03 2.5304E-7; 92235.03 6.1604E-5; @@ -84,7 +84,7 @@ nuclearData { moder {1001.03 hh2o.04; } } aluminium { - temp 293; + !temp 293; composition { 13027.03 5.9699E-2; 14000.03 5.5202E-4; diff --git a/InputFiles/Benchmarks/MCNP/ORNL10 b/InputFiles/Benchmarks/MCNP/ORNL10 index 6a201802e..d4c396646 100644 --- a/InputFiles/Benchmarks/MCNP/ORNL10 +++ b/InputFiles/Benchmarks/MCNP/ORNL10 @@ -73,7 +73,7 @@ materials { solution { - temp 293; + !temp 293; composition { 92233.03 3.9124E-9; 92234.03 4.0905E-7; @@ -87,7 +87,7 @@ materials { moder {1001.03 hh2o.04; } } shell { - temp 293; + !temp 293; composition { 13027.03 5.9881E-2; 14000.03 2.1790E-4; diff --git a/InputFiles/Benchmarks/MCNP/ORNL11 b/InputFiles/Benchmarks/MCNP/ORNL11 index f6651c976..f594aca7f 100644 --- a/InputFiles/Benchmarks/MCNP/ORNL11 +++ b/InputFiles/Benchmarks/MCNP/ORNL11 @@ -73,7 +73,7 @@ materials { solution { - temp 293; + !temp 293; composition { 92233.03 3.3441E-5; 92234.03 5.2503E-7; @@ -87,7 +87,7 @@ materials { moder { 1001.03 hh2o.04; } } shell { - temp 293; + !temp 293; composition { 13027.03 5.9881E-2; 14000.03 2.1790E-4; diff --git a/InputFiles/Benchmarks/MCNP/PNL2 b/InputFiles/Benchmarks/MCNP/PNL2 index ccdb3cea9..3378a518e 100644 --- a/InputFiles/Benchmarks/MCNP/PNL2 +++ b/InputFiles/Benchmarks/MCNP/PNL2 @@ -75,7 +75,7 @@ materials { solution { - temp 293; + !temp 293; composition { 94238.03 2.6153E-8; 94239.03 4.1249E-4; diff --git a/InputFiles/Benchmarks/MCNP/PUMF11 b/InputFiles/Benchmarks/MCNP/PUMF11 index 1d62df2d6..c1d9caade 100644 --- a/InputFiles/Benchmarks/MCNP/PUMF11 +++ b/InputFiles/Benchmarks/MCNP/PUMF11 @@ -74,7 +74,7 @@ materials { fuel { - temp 293; + !temp 293; composition { 94239.03 4.6982E-2; 94240.03 2.5852E-3; @@ -83,7 +83,7 @@ materials { } } water { - temp 293; + !temp 293; composition { 1001.03 6.6766E-2; 08016.03 3.3383E-2; diff --git a/InputFiles/Benchmarks/MCNP/SB212 b/InputFiles/Benchmarks/MCNP/SB212 index ff4c944bc..2a59c5cee 100644 --- a/InputFiles/Benchmarks/MCNP/SB212 +++ b/InputFiles/Benchmarks/MCNP/SB212 @@ -166,7 +166,7 @@ materials { fuel { // 233UO2-ZrO2 - temp 293; + !temp 293; composition { 92233.03 3.9891E-3; 92234.03 6.3690E-5; @@ -177,7 +177,7 @@ materials { } blade { - temp 293; + !temp 293; composition { 26000.03 5.9259E-2; // Fe 24000.03 1.7428E-2; // Cr @@ -188,7 +188,7 @@ materials { } water { - temp 293; + !temp 293; composition { 08016.03 3.3368E-2; 1001.03 6.6735E-2; @@ -197,7 +197,7 @@ materials { } zirc { // Zircaloy - temp 293; + !temp 293; composition { 40000.03 4.2537E-2; // Zr 50000.03 4.9918E-4; // Sn diff --git a/InputFiles/Benchmarks/MCNP/Tinkertoy2 b/InputFiles/Benchmarks/MCNP/Tinkertoy2 index 31aec3479..b79e65c73 100644 --- a/InputFiles/Benchmarks/MCNP/Tinkertoy2 +++ b/InputFiles/Benchmarks/MCNP/Tinkertoy2 @@ -140,7 +140,7 @@ nuclearData { materials { fuel { - temp 293; + !temp 293; composition { 92234.03 4.8271E-4; 92235.03 4.4797E-2; @@ -150,7 +150,7 @@ materials { } steel { - temp 293; + !temp 293; composition { 06012.03 3.1691E-4; // C 25055.03 1.7321E-3; // Mn @@ -162,7 +162,7 @@ materials { } concrete { - temp 293; + !temp 293; composition { 01001.03 1.4868E-2; // H 06012.03 3.8144E-3; // C @@ -177,7 +177,7 @@ materials { } paraffin { - temp 293; + !temp 293; composition { 01001.03 8.2574E-2; // H 06012.03 3.9699E-2; // C @@ -186,7 +186,7 @@ materials { air { // No data for air given, this composition taken from: // https://www.researchgate.net/figure/contd-Atom-Densities-for-Basic-Materials-atom-barn-cm_tbl17_267562651 - temp 293; + !temp 293; composition { 07014.03 4.1985E-5; 08016.03 1.1263E-5; diff --git a/InputFiles/Benchmarks/MCNP/U233MF05 b/InputFiles/Benchmarks/MCNP/U233MF05 index 332a1270f..25d11c59d 100644 --- a/InputFiles/Benchmarks/MCNP/U233MF05 +++ b/InputFiles/Benchmarks/MCNP/U233MF05 @@ -73,7 +73,7 @@ materials { fuel { - temp 293; + !temp 293; composition { 92233.03 0.047312; 92234.03 0.00052770; @@ -81,7 +81,7 @@ materials { } } reflector { - temp 293; + !temp 293; composition { 04009.03 0.11984; 08016.03 0.0013776; diff --git a/InputFiles/Benchmarks/SourceConvergence/checkerboard b/InputFiles/Benchmarks/SourceConvergence/checkerboard index cc01927e3..f1c79c187 100644 --- a/InputFiles/Benchmarks/SourceConvergence/checkerboard +++ b/InputFiles/Benchmarks/SourceConvergence/checkerboard @@ -186,7 +186,6 @@ nuclearData { materials { water { - temp 75675; rgb (0 0 139); // Colour of water is dark blue moder { 1001.03 lwj3.00; } composition { @@ -196,21 +195,19 @@ nuclearData { } clad { - temp 12345; composition { 40000.03 4.291E-002; } } fuel { - temp 87476; composition { 92235.03 8.2213E-004; 92238.03 2.238E-002; 8016.03 4.6054E-002; } } + concrete { - temp 6786; composition { 1001.03 5.5437E-03; 6012.03 6.9793E-03; @@ -218,8 +215,8 @@ nuclearData { 20000.03 8.9591E-03; 8016.03 4.3383E-002;} } + steel { - temp 8765; composition { 26000.03 8.377E-02; } diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 index cd197bff7..9896f9418 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 @@ -24,13 +24,13 @@ module aceNeutronDatabase_iTest ! Material definitions character(*),parameter :: MAT_INPUT_STR = & - & "water { temp 273; & + & "water { & & composition { & & 1001.03 5.028E-02; & & 8016.03 2.505E-02; & & } & & } & - & uo2 { temp 1; & + & uo2 { & & composition { & & 92233.03 2.286E-02; & & 8016.03 4.572E-02; & @@ -91,7 +91,7 @@ subroutine test_aceNeutronDatabase() @assertNotAssociated( ceNeutronMaterial_TptrCast( data % getMaterial(3))) ! Get water - mat => ceNeutronMaterial_TptrCast( data % getMaterial(1)) + mat => ceNeutronMaterial_TptrCast(data % getMaterial(1)) @assertAssociated(mat) ! Make sure densities are present @@ -172,7 +172,6 @@ subroutine test_aceNeutronDatabase() name = 'uo2' @assertTrue( 0 /= matNames % getOrDefault(name, 0)) - !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test getting nuclide XSs ! diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 index 61c75c77a..cbd3e2313 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 @@ -17,14 +17,14 @@ module thermalScatteringData_iTest ! Material definitions character(*),parameter :: MAT_INPUT_STR = & - & "water { temp 1; & + & "water { & & moder {1001.03 h-h2o.49; } & & composition { & & 1001.03 2.0E-3; & & 8016.03 1.0E-3; & & } & & } & - & graphite { temp 1; & + & graphite { & & moder {6012.06 grph30.46;} & & composition { & & 6012.06 2.0E-3; & @@ -126,9 +126,8 @@ subroutine test_thermalScatteringData() !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test getting XSs ! H-1 - nuc => ceNeutronNuclide_CptrCast( data % getNuclide(2)) + nuc => ceNeutronNuclide_CptrCast(data % getNuclide(2)) nuclideCache(2) % E_tot = ONE - nuclideCache(2) % needsSabInel = .true. call nuc % getMicroXSs(microXSs, 1.8E-6_defReal, p % pRNG) diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 index 6afb16e30..cba91b668 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 @@ -17,7 +17,7 @@ module urrProbabilityTables_iTest ! Material definitions character(*),parameter :: MAT_INPUT_STR = & - & " uo2 { temp 1; & + & " uo2 { & & composition { & & 92235.03 1.0E-3; & & 8016.03 2.0E-3; & @@ -26,7 +26,7 @@ module urrProbabilityTables_iTest ! CE Neutron Database specification character(*),parameter :: ACE_INPUT_STR = & - & "aceLibrary ./IntegrationTestFiles/testLib; ures 1 ;" + & "aceLibrary ./IntegrationTestFiles/testLib; ures 1 ; majorant 1; " contains @@ -98,10 +98,9 @@ subroutine test_urrProbabilityTables() ! U-235 nuc => ceNeutronNuclide_CptrCast( data % getNuclide(1)) - zaidCache(1) % E = 9.1E-3_defReal + zaidCache(1) % E = 9.1E-3_defReal zaidCache(1) % xi = 0.347_defReal nuclideCache(1) % E_tot = ONE - nuclideCache(1) % needsUrr = .true. call nuc % getMicroXSs(microXSs, 9.1E-3_defReal, p % pRNG) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 19190e80f..b8e36a7b2 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -37,6 +37,9 @@ module aceNeutronDatabase_class cache_zaidCache => zaidCache, & cache_init => init + ! Scattering procedures + use scatteringKernels_func, only : relativeEnergy_constXS, dopplerCorrectionFactor + implicit none private @@ -61,7 +64,7 @@ module aceNeutronDatabase_class !! Public Members: !! nuclides -> array of aceNeutronNuclides with data !! materials -> array of ceNeutronMaterials with data - !! Ebounds -> array with bottom (1) and top (2) energy bound + !! eBounds -> array with bottom (1) and top (2) energy bound !! majorant -> unionised majorant cross section !! eGridUnion -> unionised energy grid !! activeMat -> array of materials present in the geometry @@ -79,7 +82,7 @@ module aceNeutronDatabase_class type(ceNeutronMaterial),dimension(:),pointer :: materials => null() real(defReal), dimension(:), allocatable :: majorant real(defReal), dimension(:), allocatable :: eGridUnion - real(defReal), dimension(2) :: Ebounds = ZERO + real(defReal), dimension(2) :: eBounds = ZERO integer(shortInt),dimension(:),allocatable :: activeMat ! Probability tables data @@ -101,22 +104,25 @@ module aceNeutronDatabase_class ! ceNeutronDatabase Procedures procedure :: energyBounds - procedure :: updateTotalMatXS procedure :: updateMajorantXS + procedure :: updateTrackMatXS + procedure :: updateTotalMatXS procedure :: updateMacroXSs procedure :: updateTotalNucXS procedure :: updateMicroXSs + procedure :: updateTotalTempNucXS procedure :: getScattMicroMajXS ! class Procedures procedure :: initUrr procedure :: initDBRC procedure :: initMajorant + procedure :: updateTotalTempMajXS + procedure :: updateRelEnMacroXSs end type aceNeutronDatabase - contains !! @@ -136,7 +142,7 @@ elemental subroutine kill(self) deallocate(self % materials) end if - self % EBounds = ZERO + self % eBounds = ZERO if(allocated(self % activeMat)) deallocate(self % activeMat) @@ -209,37 +215,40 @@ function getReaction(self, MT, idx) result(reac) ! MT < 0 -> material reaction ! MT = 0 -> does not exist ! MT = 1 -> N_total has no reaction object - if ( MT <= 1 ) then + if (MT <= 1) then reac => null() return end if ! Detect invalid indices - if( idx < 1 .or. idx > size(self % nuclides)) then + if (idx < 1 .or. idx > size(self % nuclides)) then reac => null() return end if ! Get nuclide reaction - if( MT == N_N_elastic) then - if (cache_nuclideCache(idx) % needsSabEl) then - reac => self % nuclides(idx) % thData % elasticOut - else - reac => self % nuclides(idx) % elasticScatter - end if - else if ( MT == N_fission) then + if (MT == N_N_elastic) then + reac => self % nuclides(idx) % elasticScatter + + else if (MT == N_N_ThermEL) then + reac => self % nuclides(idx) % thData % elasticOut + + else if (MT == N_fission) then reac => self % nuclides(idx) % fission + else if (MT == N_N_ThermINEL) then reac => self % nuclides(idx) % thData % inelasticOut + else ! Find index of MT reaction idxMT = self % nuclides(idx) % idxMT % getOrDefault(MT, 0) ! See if the MT is present or not - if(idxMT == 0) then + if (idxMT == 0) then reac => null() else reac => self % nuclides(idx) % MTdata(idxMT) % kinematics end if + end if end function getReaction @@ -256,23 +265,23 @@ function getScattMicroMajXS(self, E, kT, A, nucIdx) result(maj) real(defReal), intent(in) :: A integer(shortInt), intent(in) :: nucIdx real(defReal) :: maj - real(defReal) :: E_upper, E_lower, E_min, E_max + real(defReal) :: eUpper, eLower, eMin, eMax real(defReal) :: alpha ! Find energy limits to define majorant calculation range - alpha = 4 / (sqrt( E * A / kT )) - E_upper = E * (1 + alpha) * (1 + alpha) - E_lower = E * (1 - alpha) * (1 - alpha) + alpha = 3.0_defReal * sqrt( kT / (E * A) ) + eUpper = E * (ONE + alpha) * (ONE + alpha) + eLower = E * (ONE - alpha) * (ONE - alpha) ! Find system minimum and maximum energies - call self % energyBounds(E_min, E_max) + call self % energyBounds(eMin, eMax) ! Avoid energy limits being outside system range - if (E_lower < E_min) E_lower = E_min - if (E_upper > E_max) E_upper = E_max + if (eLower < eMin .or. ONE < alpha) eLower = eMin + if (eUpper > eMax) eUpper = eMax ! Find largest elastic scattering xs in energy range given by E_lower and E_upper - maj = self % nuclides(nucIdx) % elScatteringMaj(E_lower, E_upper) + maj = self % nuclides(nucIdx) % getMajXS(eLower, eUpper, N_N_ELASTIC) end function getScattMicroMajXS @@ -281,58 +290,19 @@ end function getScattMicroMajXS !! !! See ceNeutronDatabase for more details !! - subroutine energyBounds(self, E_min, E_max) + subroutine energyBounds(self, eMin, eMax) class(aceNeutronDatabase), intent(in) :: self - real(defReal), intent(out) :: E_min - real(defReal), intent(out) :: E_max + real(defReal), intent(out) :: eMin + real(defReal), intent(out) :: eMax - E_min = self % Ebounds(1) - E_max = self % Ebounds(2) + eMin = self % eBounds(1) + eMax = self % eBounds(2) end subroutine energyBounds !! - !! Make sure that totalXS of material with matIdx is at energy E - !! in ceNeutronChache - !! - !! See ceNeutronDatabase for more details - !! - subroutine updateTotalMatXS(self, E, matIdx, rand) - class(aceNeutronDatabase), intent(in) :: self - real(defReal), intent(in) :: E - integer(shortInt), intent(in) :: matIdx - class(RNG), optional, intent(inout) :: rand - integer(shortInt) :: i, nucIdx - real(defReal) :: dens - - associate (mat => cache_materialCache(matIdx)) - ! Set new energy - mat % E_tot = E - - ! Clean current total XS - mat % xss % total = ZERO - - ! Construct total macro XS - do i = 1, size(self % materials(matIdx) % nuclides) - dens = self % materials(matIdx) % dens(i) - nucIdx = self % materials(matIdx) % nuclides(i) - - ! Update if needed - if (cache_nuclideCache(nucIdx) % E_tot /= E) then - call self % updateTotalNucXS(E, nucIdx, rand) - end if - - ! Add microscopic XSs - mat % xss % total = mat % xss % total + dens * cache_nuclideCache(nucIdx) % xss % total - end do - - end associate - - end subroutine updateTotalMatXS - - !! - !! Make sure that the majorant of ALL Active materials is at energy E - !! in ceNeutronChache + !! Make sure that the majorant of ALL active materials is at energy E + !! in ceNeutronCache !! !! See ceNeutronDatabase for more details !! @@ -344,20 +314,20 @@ subroutine updateMajorantXS(self, E, rand) real(defReal) :: f character(100), parameter :: Here = 'updateMajorantXS (aceNeutronDatabase_class.f90)' - associate (maj => cache_majorantCache(1) ) + associate (maj => cache_majorantCache(1)) maj % E = E ! Get majorant via the precomputed unionised cross section if (self % hasMajorant) then idx = binarySearch(self % eGridUnion, E) - if(idx <= 0) then + if (idx <= 0) then call fatalError(Here,'Failed to find energy: '//numToChar(E)//& ' in unionised majorant grid') end if - associate(E_top => self % eGridUnion(idx + 1), E_low => self % eGridUnion(idx)) + associate (E_top => self % eGridUnion(idx + 1), E_low => self % eGridUnion(idx)) f = (E - E_low) / (E_top - E_low) end associate @@ -372,11 +342,11 @@ subroutine updateMajorantXS(self, E, rand) matIdx = self % activeMat(i) ! Update if needed - if( cache_materialCache(matIdx) % E_tot /= E) then - call self % updateTotalMatXS(E, matIdx, rand) + if (cache_materialCache(matIdx) % E_track /= E) then + call self % updateTrackMatXS(E, matIdx, rand) end if - maj % xs = max(maj % xs, cache_materialCache(matIdx) % xss % total) + maj % xs = max(maj % xs, cache_materialCache(matIdx) % trackXS) end do end if @@ -385,6 +355,132 @@ subroutine updateMajorantXS(self, E, rand) end subroutine updateMajorantXS + !! + !! Make sure that trackXS of material with matIdx is at energy E = E_track + !! in ceNeutronChache + !! + !! See ceNeutronDatabase for more details + !! + subroutine updateTrackMatXS(self, E, matIdx, rand) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + class(RNG), optional, intent(inout) :: rand + + associate (matCache => cache_materialCache(matIdx), & + mat => self % materials(matIdx)) + + ! Set new energy + matCache % E_track = E + + if (mat % useTMS(E)) then + ! The material tracking xs is the temperature majorant in the case of TMS + call self % updateTotalTempMajXS(E, matIdx) + + else + ! When TMS is not in use, the material tracking xs is equivalent to the total + call self % updateTotalMatXS(E, matIdx, rand) + matCache % trackXS = matCache % xss % total + + end if + + end associate + + end subroutine updateTrackMatXS + + !! + !! Subroutine to update the temperature majorant in a given material at given temperature + !! + !! The function finds the upper and lower limits of the energy range the nuclide majorant + !! is included into, then adds up the nuclide temperature majorant multiplied by the Doppler + !! correction factor + !! + !! Args: + !! E [in] -> Incident neutron energy for which temperature majorant is found + !! matIdx [in] -> Index of material for which the material temperature majorant is found + !! + subroutine updateTotalTempMajXS(self, E, matIdx) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + integer(shortInt) :: nucIdx, i + real(defReal) :: dens, corrFact, nucTempMaj + + associate (matCache => cache_materialCache(matIdx), & + mat => self % materials(matIdx)) + + ! Clean current total XS + matCache % trackXS = ZERO + + ! loop through all nuclides in material and find sum of majorants + do i = 1, size(mat % nuclides) + + ! Get nuclide data + nucIdx = mat % nuclides(i) + dens = mat % dens(i) + + call self % updateTotalTempNucXS(E, mat % kT, nucIdx) + + ! Sum nuclide majorants to find material majorant + corrFact = cache_nuclideCache(nucIdx) % doppCorr + nucTempMaj = cache_nuclideCache(nucIdx) % tempMajXS * corrFact + matCache % trackXS = matCache % trackXS + dens * nucTempMaj + + end do + + end associate + + end subroutine updateTotalTempMajXS + + !! + !! Make sure that totalXS of material with matIdx is at energy E + !! in ceNeutronCache + !! + !! See ceNeutronDatabase for more details + !! + subroutine updateTotalMatXS(self, E, matIdx, rand) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + class(RNG), optional, intent(inout) :: rand + integer(shortInt) :: i, nucIdx + real(defReal) :: dens + + associate (matCache => cache_materialCache(matIdx), & + mat => self % materials(matIdx)) + + ! Set new energy and clean current total XS + matCache % E_tot = E + matCache % xss % total = ZERO + + if (mat % useTMS(E)) then + ! When TMS is in use, the total xs is retrieved sampling the nuclides' relative + ! energies given the temperature difference between material temperature and + ! temperature of the nuclides' base cross sections + call self % updateRelEnMacroXSs(E, matIdx, rand) + + else + ! Construct total macro XS + do i = 1, size(mat % nuclides) + dens = mat % dens(i) + nucIdx = mat % nuclides(i) + + ! Update if needed + if (cache_nuclideCache(nucIdx) % E_tot /= E) then + call self % updateTotalNucXS(E, nucIdx, rand) + end if + + ! Add microscopic XSs + matCache % xss % total = matCache % xss % total + & + dens * cache_nuclideCache(nucIdx) % xss % total + end do + + end if + + end associate + + end subroutine updateTotalMatXS + !! !! Make sure that the macroscopic XSs for the material with matIdx are set !! to energy E in ceNeutronCache @@ -399,35 +495,122 @@ subroutine updateMacroXSs(self, E, matIdx, rand) integer(shortInt) :: i, nucIdx real(defReal) :: dens - associate (mat => cache_materialCache(matIdx)) - ! Set new energy - mat % E_tot = E - mat % E_tail = E + associate(mat => self % materials(matIdx), & + matCache => cache_materialCache(matIdx)) ! Clean current xss - call mat % xss % clean() + call matCache % xss % clean() - ! Construct microscopic XSs - do i = 1, size(self % materials(matIdx) % nuclides) - dens = self % materials(matIdx) % dens(i) - nucIdx = self % materials(matIdx) % nuclides(i) + if (mat % useTMS(E)) then + ! When TMS is in use, the xss are retrieved sampling the nuclides' relative + ! energies given the temperature difference between material temperature and + ! temperature of the nuclides' base cross sections + call self % updateRelEnMacroXSs(E, matIdx, rand) - ! Update if needed - if (cache_nuclideCache(nucIdx) % E_tail /= E .or. cache_nuclideCache(nucIdx) % E_tot /= E) then - call self % updateMicroXSs(E, nucIdx, rand) - end if + else - ! Add microscopic XSs - call mat % xss % add(cache_nuclideCache(nucIdx) % xss, dens) - end do + ! Set new energy + matCache % E_tot = E + matCache % E_tail = E + + ! Construct microscopic XSs + do i = 1, size(mat % nuclides) + dens = mat % dens(i) + nucIdx = mat % nuclides(i) + + ! Update if needed + if (cache_nuclideCache(nucIdx) % E_tail /= E .or. cache_nuclideCache(nucIdx) % E_tot /= E) then + call self % updateMicroXSs(E, nucIdx, rand) + end if + + ! Add microscopic XSs + call matCache % xss % add(cache_nuclideCache(nucIdx) % xss, dens) + end do + + end if end associate end subroutine updateMacroXSs + !! + !! Subroutine to update the macroscopic cross sections in a given material + !! at given temperature, sampling a relative energy per each nuclide and applying + !! a Doppler correction factor + !! + !! Args: + !! E [in] -> Incident neutron energy for which the relative energy xss are found + !! matIdx [in] -> Index of material for which the relative energy xss are found + !! + subroutine updateRelEnMacroXSs(self, E, matIdx, rand) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + class(RNG), optional, intent(inout) :: rand + integer(shortInt) :: i, nucIdx + real(defReal) :: dens, nuckT, A, deltakT, eRel, eMin, & + eMax, doppCorr + character(100), parameter :: Here = 'updateRelEnMacroXSs (aceNeutronDatabase_class.f90)' + + associate(mat => self % materials(matIdx), & + matCache => cache_materialCache(matIdx)) + + ! Check if relative energy cross sections have been retrieved before + if (E /= matCache % E_rel) then + + ! Clean current xss + call matCache % xssRel % clean() + matCache % E_rel = E + + ! Construct microscopic XSs + do i = 1, size(mat % nuclides) + + dens = mat % dens(i) + nucIdx = mat % nuclides(i) + nuckT = self % nuclides(nucIdx) % getkT() + A = self % nuclides(nucIdx) % getMass() + deltakT = mat % kT - nuckT + + eRel = relativeEnergy_constXS(E, A, deltakT, rand) + + ! Call through system minimum and maximum energies + call self % energyBounds(eMin, eMax) + + ! avoid sampled relative energy from MB dist extending into energies outside system range + if (eRel < eMin) eRel = eMin + if (eMax < eRel) eRel = eMax + + associate(nucCache => cache_nuclideCache(nucIdx)) + + ! Doppler correction factor for low energies + doppCorr = dopplerCorrectionFactor(E, A, deltakT) + + ! Update if needed + if (nucCache % E_tail /= eRel .or. nucCache % E_tot /= eRel) then + call self % updateMicroXSs(eRel, nucIdx, rand) + end if + + ! Add microscopic XSs + call matCache % xssRel % add(nucCache % xss, dens * doppCorr) + + end associate + + end do + + end if + + ! Update cache, and ensure that the energy indicators are reset to avoid wrong look-ups + matCache % xss = matCache % xssRel + matCache % E_tot = ZERO + matCache % E_tail = ZERO + + end associate + + end subroutine updateRelEnMacroXSs + !! !! Make sure that totalXS of nuclide with nucIdx is at energy E - !! in ceNeutronChache + !! in ceNeutronCache !! !! See ceNeutronDatabase for more details !! @@ -436,19 +619,12 @@ subroutine updateTotalNucXS(self, E, nucIdx, rand) real(defReal), intent(in) :: E integer(shortInt), intent(in) :: nucIdx class(RNG), optional, intent(inout) :: rand - logical(defBool) :: needsSab associate (nucCache => cache_nuclideCache(nucIdx), & nuc => self % nuclides(nucIdx) ) - ! Check if the nuclide needs ures probability tables at this energy - nucCache % needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2)) - ! Check if the nuclide needs S(a,b) at this energy - nucCache % needsSabEl = (nuc % hasThData .and. E >= nuc % SabEl(1) .and. E <= nuc % SabEl(2)) - nucCache % needsSabInel = (nuc % hasThData .and. E >= nuc % SabInel(1) .and. E <= nuc % SabInel(2)) - needsSab = (nucCache % needsSabEl .or. nucCache % needsSabInel) - - if (nucCache % needsUrr .or. needsSab) then + ! Check if the nuclide needs ures probability tables or S(a,b) at this energy + if (nuc % needsUrr(E) .or. nuc % needsSabEl(E) .or. nuc % needsSabInel(E)) then call self % updateMicroXSs(E, nucIdx, rand) else @@ -487,7 +663,7 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) ! Overwrites all the micro cross sections in cache ! Check if probability tables should be read - if (nucCache % needsUrr) then + if (nuc % needsUrr(E)) then associate(zaidCache => cache_zaidCache(self % nucToZaid(nucIdx))) if (zaidCache % E /= E) then @@ -500,7 +676,8 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) end associate - elseif (nucCache % needsSabEl .or. nucCache % needsSabInel) then + ! Check if S(a,b) should be read + elseif (nuc % needsSabEl(E) .or. nuc % needsSabInel(E)) then call nuc % getThXSs(nucCache % xss, nucCache % idx, nucCache % f, E) else @@ -512,6 +689,60 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) end subroutine updateMicroXSs + !! + !! Subroutine to retrieve the nuclide total majorant cross section + !! over a calculated energy range (needed for TMS) and update the nuclide + !! cache with majorant value, energy, deltakT and Doppler correction factor + !! + !! See ceNeutronDatabase for more details + !! + subroutine updateTotalTempNucXS(self, E, kT, nucIdx) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + real(defReal), intent(in) :: kT + integer(shortInt), intent(in) :: nucIdx + real(defReal) :: eUpper, eLower, eMin, eMax, nuckT, & + alpha, deltakT, A + character(100), parameter :: Here = 'updateTotalTempNucXS (aceNeutronDatabase_class.f90)' + + associate (nuc => self % nuclides(nucIdx) , & + nucCache => cache_nuclideCache(nucIdx)) + + nuckT = nuc % getkT() + A = nuc % getMass() + deltakT = kT - nuckT + + ! Check if an update is required + if (nucCache % E_maj /= E .or. nucCache % deltakT /= deltakT) then + + ! Find energy limits to define majorant calculation range + alpha = 3.0_defReal * sqrt( deltakT / (E * A) ) + eUpper = E * (ONE + alpha) * (ONE + alpha) + eLower = E * (ONE - alpha) * (ONE - alpha) + + ! Find system minimum and maximum energies + call self % energyBounds(eMin, eMax) + + ! Avoid energy limits being outside system range + if (eLower < eMin .or. ONE < alpha) eLower = eMin + if (eUpper > eMax) eUpper = eMax + + ! Doppler g correction factor for low energies + nucCache % doppCorr = dopplerCorrectionFactor(E, A, deltakT) + + ! Get nuclide total majorant cross section in the energy range + nucCache % tempMajXS = nuc % getMajXS(eLower, eUpper, N_TOTAL) + + ! Save additional info + nucCache % deltakT = deltakT + nucCache % E_maj = E + + end if + + end associate + + end subroutine updateTotalTempNucXS + !! !! Initialise Database from dictionary and pointer to self !! @@ -536,6 +767,9 @@ subroutine init(self, dict, ptr, silent ) logical(defBool) :: isFissileMat integer(shortInt),dimension(:),allocatable :: nucIdxs, zaidDBRC character(nameLen),dimension(:),allocatable :: nucDBRC + real(defReal) :: A, nuckT, eUpSab, eUpSabNuc, & + eLowURR, eLowUrrNuc, alpha, & + deltakT, eUpper, eLower integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 character(100), parameter :: Here = 'init (aceNeutronDatabase_class.f90)' @@ -660,6 +894,17 @@ subroutine init(self, dict, ptr, silent ) i = nucSet % next(i) end do + ! Calculate energy bounds + self % eBounds(1) = self % nuclides(1) % eGrid(1) + j = size(self % nuclides(1) % eGrid) + self % eBounds(2) = self % nuclides(1) % eGrid(j) + + do i = 2, size(self % nuclides) + self % eBounds(1) = max(self % eBounds(1), self % nuclides(i) % eGrid(1)) + j = size(self % nuclides(i) % eGrid) + self % eBounds(2) = min(self % eBounds(2), self % nuclides(i) % eGrid(j)) + end do + ! Build Material definitions allocate(self % materials(mm_nMat())) allocate(nucIdxs(maxNuc)) @@ -667,35 +912,82 @@ subroutine init(self, dict, ptr, silent ) mat => mm_getMatPtr(i) ! Load nuclide indices on storage space + ! Find if material is fissile isFissileMat = .false. + ! Loop over nuclides do j = 1, size(mat % nuclides) name = trim(mat % nuclides(j) % toChar()) + if (mat % nuclides(j) % hasSab) then zaid = trim(mat % nuclides(j) % toChar()) file = trim(mat % nuclides(j) % file_Sab) name = zaid // '+' // file deallocate(zaid, file) end if + + ! Find nuclide definition to see if fissile nucIdxs(j) = nucSet % get(name) isFissileMat = isFissileMat .or. self % nuclides(nucIdxs(j)) % isFissile() + end do ! Load data into material - call self % materials(i) % set( matIdx = i, & + call self % materials(i) % set( matIdx = i, & database = ptr_ceDatabase, & - fissile = isFissileMat ) + temp = mat % T, & + hasTMS = mat % hasTMS, & + fissile = isFissileMat ) call self % materials(i) % setComposition( mat % dens, nucIdxs(1:size(mat % nuclides))) - end do - ! Calculate energy bounds - self % Ebounds(1) = self % nuclides(1) % eGrid(1) - j = size(self % nuclides(1) % eGrid) - self % Ebounds(2) = self % nuclides(1) % eGrid(j) + eUpSab = self % eBounds(1) + eLowURR = self % eBounds(2) + + if (mat % hasTMS) then + + ! Loop again to find energy limits of S(a,b) and URES for TMS applicability + do j = 1, size(mat % nuclides) + + ! Find nuclide information + idx = nucIdxs(j) + nuckT = self % nuclides(idx) % getkT() + A = self % nuclides(idx) % getMass() + deltakT = self % materials(i) % kT - nuckT + + ! Call fatal error if material temperature is lower then base nuclide temperature + if (deltakT < ZERO) then + call fatalError(Here, "Material temperature must be greater than the nuclear data temperature.") + end if + + ! Find nuclide upper S(a,b) energy + eUpSabNuc = max(self % nuclides(idx) % SabEl(2), self % nuclides(idx) % SabInel(2)) + + ! Find energy limits to define majorant calculation range + if (eUpSabNuc > ZERO) then + alpha = 4.0_defReal * sqrt( deltakT / (eUpSabNuc * A) ) + eUpper = eUpSabNuc * (ONE + alpha) * (ONE + alpha) + else + eUpper = ZERO + end if + + eLowUrrNuc = self % nuclides(idx) % urrE(1) + + if (eLowUrrNuc /= ZERO) then + alpha = 4.0_defReal * sqrt( deltakT / (eLowUrrNuc * A) ) + eLower = eLowUrrNuc * (ONE - alpha) * (ONE - alpha) + else + eLower = self % eBounds(2) + end if + + eUpSab = max(eUpSab, eUpper) + eLowURR = min(eLowURR, eLower) + + end do + + end if + + ! Load data into material + call self % materials(i) % set(eUpperSab = eUpSab, eLowerURR = eLowURR) - do i = 2, size(self % nuclides) - self % Ebounds(1) = max(self % Ebounds(1), self % nuclides(i) % eGrid(1)) - j = size(self % nuclides(i) % eGrid) - self % Ebounds(2) = min(self % Ebounds(2), self % nuclides(i) % eGrid(j)) end do ! Read unionised majorant flag @@ -821,7 +1113,7 @@ subroutine activate(self, activeMat, silent) ! Configure Cache if (self % hasUrr) then - call cache_init(size(self % materials), size(self % nuclides), 1, maxval(self % nucToZaid)) + call cache_init(size(self % materials), size(self % nuclides), nZaid = maxval(self % nucToZaid)) else call cache_init(size(self % materials), size(self % nuclides)) end if @@ -856,9 +1148,9 @@ subroutine initMajorant(self, loud) sizeGrid, eIdx, nucIdxLast, eIdxLast, & urrIdx type(intMap) :: nucSet - real(defReal) :: eRef, eNuc, E, maj, total, dens, urrMaj, & + real(defReal) :: eRef, eNuc, E, maj, trackXS, dens, urrMaj, & nucXS, f, eMax, eMin - logical(defBool) :: needsUrr + class(RNG), allocatable :: rand integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 real(defReal), parameter :: NUDGE = 1.0e-06_defReal @@ -901,7 +1193,7 @@ subroutine initMajorant(self, loud) ! Allocate temporary grid vector and initialise to largest value allowed allocate(tmpGrid(sizeGrid)) - tmpGrid = self % EBounds(2) + tmpGrid = self % eBounds(2) ! Loop over the energy grid i = 1 @@ -926,7 +1218,7 @@ subroutine initMajorant(self, loud) eNuc = self % nuclides(nucIdx) % eGrid(eIdx) ! Check if the energy from the nuclide grid is out of bounds - if (eNuc < self % EBounds(1) .or. eNuc > self % EBounds(2)) then + if (eNuc < self % eBounds(1) .or. eNuc > self % eBounds(2)) then j = nucSet % next(j) cycle end if @@ -1031,6 +1323,12 @@ subroutine initMajorant(self, loud) ! Allocate unionised majorant allocate(self % majorant(size(self % eGridUnion))) + ! Initialise RNG needed to call update XS routines. The initial seed doesn't + ! matter because the RNG is only used to sample from probability tables, which + ! are corrected for the majorant anyway + allocate(rand) + call rand % init(1_longInt) + ! Loop over all the energies do i = 1, size(self % eGridUnion) @@ -1038,8 +1336,8 @@ subroutine initMajorant(self, loud) E = self % eGridUnion(i) ! Correct for energies higher or lower than the allowed boundaries - if (E < self % EBounds(1)) E = self % EBounds(1) - if (E > self % EBounds(2)) E = self % EBounds(2) + if (E < self % eBounds(1)) E = self % eBounds(1) + if (E > self % eBounds(2)) E = self % eBounds(2) ! Initialise majorant value for this energy maj = ZERO @@ -1049,19 +1347,20 @@ subroutine initMajorant(self, loud) ! Get material index matIdx = self % activeMat(j) - total = ZERO - ! Loop over nuclides + ! Get material tracking cross section + call self % updateTrackMatXS(E, matIdx, rand) + trackXS = cache_materialCache(matIdx) % trackXS + + ! Loop over nuclides to check and correct for ures do k = 1, size(self % materials(matIdx) % nuclides) dens = self % materials(matIdx) % dens(k) nucIdx = self % materials(matIdx) % nuclides(k) associate (nuc => self % nuclides(nucIdx)) - needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2)) - ! Check if present nuclide uses URR tables - if (needsUrr) then + if (nuc % needsUrr(E)) then ! Find maximum URR table total XS urrIdx = binarySearch(nuc % probTab % eGrid, E) @@ -1070,25 +1369,22 @@ subroutine initMajorant(self, loud) ! Check if URR tables contain xs or multiplicative factor if (nuc % IFF == 1) then call nuc % search(eIdx, f, E) - nucXS = nuc % totalXS(eIdx, f) * urrMaj + nucXS = nuc % totalXS(eIdx, f) * urrMaj else nucXS = urrMaj end if - else - call self % updateTotalNucXS(E, nucIdx) - nucXS = cache_nuclideCache(nucIdx) % xss % total + ! Update total material cross section + trackXS = trackXS + dens * (nucXS - cache_nuclideCache(nucIdx) % xss % total) end if end associate - ! Update total material cross section - total = total + dens * nucXS - end do - maj = max(maj, total) + ! Select majorant cross section + maj = max(maj, trackXS) end do diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 343bc6b0e..5715e645c 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -101,7 +101,7 @@ module aceNeutronNuclide_class !! microXSs -> return interpolated ceNeutronMicroXSs package given index and inter. factor !! getUrrXSs -> return ceNeutronMicroXSs accounting for ures probability tables !! getThXSs -> return ceNeutronMicroXSs accounting for S(a,b) scattering treatment - !! elScatteringMaj -> returns the elastic scattering majorant within an energy range given as input + !! getMajXS -> returns a majorant cross section on request within an energy range given as input !! init -> build nuclide from aceCard !! initUrr -> build list and mapping of nuclides to maintain temperature correlation !! when reading ures probability tables @@ -136,6 +136,7 @@ module aceNeutronNuclide_class procedure :: invertInelastic procedure :: xsOf procedure :: elScatteringXS + procedure :: needsSabEl procedure :: kill ! Local interface @@ -145,7 +146,9 @@ module aceNeutronNuclide_class procedure :: microXSs procedure :: getUrrXSs procedure :: getThXSs - procedure :: elScatteringMaj + procedure :: getMajXS + procedure :: needsUrr + procedure :: needsSabInel procedure :: init procedure :: initUrr procedure :: initSab @@ -174,29 +177,24 @@ function invertInelastic(self, E, rand) result(MT) character(100), parameter :: Here = 'invertInelastic (aceNeutronNuclide_class.f90)' ! Check if it's thermal inelastic scattering or not - if (nuclideCache(self % getNucIdx()) % needsSabInel) then + if (self % needsSabInel(E)) then MT = N_N_ThermINEL return end if ! Normal (without S(a,b)) inelastic scattering ! Obtain bin index and interpolation factor - if (nuclideCache(self % getNucIdx()) % E_tot == E) then - idx = nuclideCache(self % getNucIdx()) % idx - f = nuclideCache(self % getNucIdx()) % f - else - call self % search(idx, f, E) - end if + call self % search(idx, f, E) ! Get inelastic XS XS = self % mainData(IESCATTER_XS, idx+1) * f + (ONE-f) * self % mainData(IESCATTER_XS, idx) ! Invert XS = XS * rand % get() - do i=1,self % nMT + do i = 1,self % nMT ! Get index in MT reaction grid idxT = idx - self % MTdata(i) % firstIdx + 1 - if( idxT < 1 ) cycle + if ( idxT < 1 ) cycle ! Get top and bottom XS topXS = self % MTdata(i) % xs(idxT+1) @@ -204,7 +202,7 @@ function invertInelastic(self, E, rand) result(MT) ! Decrement total inelastic and exit if sampling is finished XS = XS - topXS * f - (ONE-f) * bottomXS - if(XS <= ZERO) then + if (XS <= ZERO) then MT = self % MTdata(i) % MT return end if @@ -489,7 +487,7 @@ elemental subroutine getThXSs(self, xss, idx, f, E) xss % elasticScatter = self % thData % getElXS(E) ! If ineleastic scatter is on, reads S(a,b) tables for inelastic scatter - if (nuclideCache(self % getNucIdx()) % needsSabInel) then + if (self % needsSabInel(E)) then xss % inelasticScatter = self % thData % getInelXS(E) else xss % inelasticScatter = data(IESCATTER_XS, 2) * f + (ONE-f) * data(IESCATTER_XS, 1) @@ -589,24 +587,40 @@ end subroutine getUrrXSs !! an energy range given by an upper and lower energy bound. !! !! Args: - !! upperE [in] -> Upper bound of energy range - !! upperE [in] -> Upper bound of energy range + !! eLower [in] -> Lower bound of energy range + !! eUpper [in] -> Upper bound of energy range + !! MT [in] -> MT number of the requested reaction cross section !! maj [out] -> Maximum scattering cross section within energy range !! - function elScatteringMaj(self, lowerE, upperE) result (maj) + function getMajXS(self, eLower, eUpper, MT) result (maj) class(aceNeutronNuclide), intent(in) :: self - real(defReal), intent(in) :: lowerE - real(defReal), intent(in) :: upperE + real(defReal), intent(in) :: eLower + real(defReal), intent(in) :: eUpper + integer(shortInt), intent(in) :: MT real(defReal) :: maj - integer(shortInt) :: idx + integer(shortInt) :: reaction, idx real(defReal) :: f, E, xs + character(100), parameter :: Here = 'getMajXS (aceNeutronNuclide_class.f90)' + + ! Select desired reaction based on requested MT number + select case (MT) + + case (N_TOTAL) + reaction = TOTAL_XS + + case (N_N_ELASTIC) + reaction = ESCATTER_XS + + case default + call fatalError(Here, 'Unsupported MT number requested: '//numToChar(MT)) + + end select ! Search for idx, f, and xs for the lower energy limit - call self % search(idx, f, lowerE) + call self % search(idx, f, eLower) ! Conservative: choose the xs at the energy point before the lower energy limit - f = 0 - maj = self % scatterXS(idx, f) + maj = self % mainData(reaction, idx) majorantLoop: do @@ -614,20 +628,69 @@ function elScatteringMaj(self, lowerE, upperE) result (maj) idx = idx + 1 ! Find XS and energy at index - xs = self % mainData(ESCATTER_XS, idx) - E = self % eGrid(idx) + xs = self % mainData(reaction, idx) + E = self % eGrid(idx) ! Compare cross sections and possibly update majorant - if (xs > maj) then - maj = xs - end if + maj = max(xs, maj) ! Exit loop after getting to the upper energy limit - if (E > upperE) exit majorantLoop + if (E >= eUpper) exit majorantLoop end do majorantLoop - end function elScatteringMaj + end function getMajXS + + !! + !! Function that checks whether this nuclide at the provided energy should + !! read unresolved resonance probability tables or not + !! + !! Args: + !! E [in] -> incident neutron energy + !! + !! Returns true or false + !! + elemental function needsUrr(self, E) result(doesIt) + class(aceNeutronNuclide), intent(in) :: self + real(defReal), intent(in) :: E + logical(defBool) :: doesIt + + doesIt = self % hasProbTab .and. E >= self % urrE(1) .and. E <= self % urrE(2) + + end function needsUrr + + !! + !! Function that checks whether this nuclide at the provided energy should + !! has S(a,b) inelastic scattering data or not + !! + !! Args: + !! E [in] -> incident neutron energy + !! + !! Returns true or false + !! + elemental function needsSabInel(self, E) result(doesIt) + class(aceNeutronNuclide), intent(in) :: self + real(defReal), intent(in) :: E + logical(defBool) :: doesIt + + doesIt = self % hasThData .and. E >= self % SabInel(1) .and. E <= self % SabInel(2) + + end function needsSabInel + + !! + !! Function that checks whether this nuclide at the provided energy should + !! has S(a,b) elastic scattering data or not + !! + !! See ceNeutronNuclide documentation + !! + elemental function needsSabEl(self, E) result(doesIt) + class(aceNeutronNuclide), intent(in) :: self + real(defReal), intent(in) :: E + logical(defBool) :: doesIt + + doesIt = self % hasThData .and. E >= self % SabEl(1) .and. E <= self % SabEl(2) + + end function needsSabEl !! !! Initialise from an ACE Card @@ -785,11 +848,11 @@ subroutine init(self, ACE, nucIdx, database) ! Calculate Inelastic scattering XS do i = 1,self % nMT - do j=1,size(self % mainData, 2) + do j = 1,size(self % mainData, 2) ! Find bottom and Top of the grid bottom = self % MTdata(i) % firstIdx top = size(self % MTdata(i) % xs) - if( j>= bottom .and. j <= top + bottom) then + if (j>= bottom .and. j <= top + bottom) then self % mainData(IESCATTER_XS, j) = self % mainData(IESCATTER_XS, j) + & self % MTdata(i) % xs(j-bottom + 1) end if @@ -797,7 +860,7 @@ subroutine init(self, ACE, nucIdx, database) end do ! Recalculate totalXS - if(self % isFissile()) then + if (self % isFissile()) then K = FISSION_XS else K = CAPTURE_XS @@ -814,7 +877,7 @@ subroutine init(self, ACE, nucIdx, database) call self % idxMT % add(N_N_ELASTIC, -ESCATTER_XS) call self % idxMT % add(N_DISAP, -CAPTURE_XS) - if(self % isFissile()) then + if (self % isFissile()) then call self % idxMT % add(N_FISSION, -FISSION_XS) end if diff --git a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 index 345b96b39..2806f3cf2 100644 --- a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 +++ b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 @@ -33,7 +33,11 @@ module ceNeutronCache_mod !! E_tail -> Energy of all XSs in xss except total !! f -> Interpolation factor for the nuclide at energy E_tot !! idx -> Index on a nuclide grid for energy E_tot - !! xss -> Cached Cross-Section values + !! xss -> Cached cross-section values + !! E_track -> Energy of the tracking xs + !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS + !! E_rel -> Base energy for which relative energy cross sections are found (for TMS) + !! xssRel -> Cached effective cross-section values at energy relative to E_rel (for TMS) !! type, public :: cacheMatDat real(defReal) :: E_tot = ZERO @@ -41,6 +45,15 @@ module ceNeutronCache_mod real(defReal) :: f = ZERO integer(shortInt) :: idx = 0 type(neutronMacroXSs) :: xss + + ! Tracking data + real(defReal) :: E_track = ZERO + real(defReal) :: trackXS = ZERO + + ! TMS data + real(defReal) :: E_rel = ZERO + type(neutronMacroXSs) :: xssRel + end type cacheMatDat !! @@ -52,12 +65,10 @@ module ceNeutronCache_mod !! f -> Interpolation factor for the nuclide at energy E_tot !! idx -> Index on a nuclide grid for energy E_tot !! xss -> Cached Cross-Sections values - !! needsSabInel -> Flag that tells if the nuclide is using thermal inelastic - !! scattering data - !! needsSabEl -> Flag that tells if the nuclide is using thermal elastic - !! scattering data - !! needsUrr -> Flag that tells if the nuclide is using unresolved resonance - !! probability tables + !! E_maj -> Energy at which the nuclide temperature majorant xs is stored (for TMS) + !! deltakT -> Difference between TMS material and nuclide thermal energy (for TMS) [MeV] + !! tempMajXS -> Temperature majorant xs value (for TMS) + !! doppCorr -> Doppler correction factor value (for TMS) !! type, public :: cacheNucDat real(defReal) :: E_tot = ZERO @@ -65,9 +76,13 @@ module ceNeutronCache_mod real(defReal) :: f = ZERO integer(shortInt) :: idx = 0 type(neutronMicroXSs) :: xss - logical(defBool) :: needsSabInel = .false. - logical(defBool) :: needsSabEl = .false. - logical(defBool) :: needsUrr = .false. + + ! TMS data + real(defReal) :: E_maj = ZERO + real(defReal) :: deltakT = ZERO + real(defReal) :: tempMajXS = ZERO + real(defReal) :: doppCorr = ONE + end type cacheNucDat !! @@ -78,8 +93,8 @@ module ceNeutronCache_mod !! xs -> value of the cross section !! type, public :: cacheSingleXS - real(defReal) :: E - real(defReal) :: xs + real(defReal) :: E = ZERO + real(defReal) :: xs = ZERO end type cacheSingleXS !! @@ -140,10 +155,10 @@ subroutine init(nMat, nNuc, nMaj, nZaid) nLoc = 1 end if - ! Chack the provided data + ! Check the provided data if (nMat < 1) call fatalError(Here,'Number of materials must be +ve! Not: '//numToChar(nMat)) - if (nNuc < 1) call fatalError(Here,'Number of nuclides must be +ve! Not: '//numToChar(nMat)) - if (nLoc < 1) call fatalError(Here,'Number of majorant XSs must be +ve! Not: '//numToChar(nMat)) + if (nNuc < 1) call fatalError(Here,'Number of nuclides must be +ve! Not: '//numToChar(nNuc)) + if (nLoc < 1) call fatalError(Here,'Number of majorant XSs must be +ve! Not: '//numToChar(nLoc)) ! Allocate space ! Need to do in parallel region to allocate each copy @@ -168,6 +183,7 @@ end subroutine init !! Return Cache Module (Singleton) to uninitialised state !! subroutine kill() + ! Need to deallocate on all threads !$omp parallel if (allocated(materialCache)) deallocate (materialCache) @@ -176,6 +192,7 @@ subroutine kill() if (allocated(trackingCache)) deallocate (trackingCache) if (allocated(zaidCache)) deallocate (zaidCache) !$omp end parallel + end subroutine kill diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index 1050c5d5b..2f9eb99b9 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -45,48 +45,75 @@ module ceNeutronDatabase_inter !! updateMacroXSs -> update Macroscopic XSs for a selected material !! updateTotalXS -> update Total XS for a selected nuclide !! updateMicroXSs -> update Microscopic XSs for a selected nuclide - !! getScattMicroMajXS -> returns elastic scattering microscopic xs majorant + !! getScattMicroMajXS -> returns nuclide elastic scattering temperature majorant (for DBRC) + !! updateTotalTempNucXS -> returns nuclide total temperature majorant (for TMS) !! type, public, abstract, extends(nuclearDatabase) :: ceNeutronDatabase type(intMap) :: mapDBRCnuc contains + ! nuclearDatabase Interface Implementation procedure :: getTrackingXS + procedure :: getTrackMatXS procedure :: getTotalMatXS procedure :: getMajorantXS ! Procedures implemented by a specific CE Neutron Database - procedure(updateTotalMatXS),deferred :: updateTotalMatXS - procedure(updateMajorantXS),deferred :: updateMajorantXS - procedure(updateMacroXSs),deferred :: updateMacroXSs - procedure(updateTotalXS),deferred :: updateTotalNucXS - procedure(updateMicroXSs),deferred :: updateMicroXSs - procedure(energyBounds),deferred :: energyBounds - procedure(getScattMicroMajXS),deferred :: getScattMicroMajXS + procedure(updateMajorantXS), deferred :: updateMajorantXS + procedure(updateTotalMatXS), deferred :: updateTrackMatXS + procedure(updateTotalMatXS), deferred :: updateTotalMatXS + procedure(updateMacroXSs), deferred :: updateMacroXSs + procedure(updateTotalXS), deferred :: updateTotalNucXS + procedure(updateMicroXSs), deferred :: updateMicroXSs + procedure(updateTotalTempNucXS), deferred :: updateTotalTempNucXS + procedure(energyBounds), deferred :: energyBounds + procedure(getScattMicroMajXS), deferred :: getScattMicroMajXS end type ceNeutronDatabase abstract interface !! !! Return energy bounds for data in the database !! - !! E_min and E_max are minimun and maximumum energy such that data + !! eMin and eMax are minimun and maximumum energy such that data !! for ALL nuclides if avalible !! !! Args: - !! E_min [out] -> minimum value of energy [MeV] - !! E_max [out] -> maximum value of energy [MeV] + !! eMin [out] -> minimum value of energy [MeV] + !! eMax [out] -> maximum value of energy [MeV] !! !! Errors: !! None !! - subroutine energyBounds(self, E_min, E_max) + subroutine energyBounds(self, eMin, eMax) import :: ceNeutronDatabase, defReal class(ceNeutronDatabase), intent(in) :: self - real(defReal), intent(out) :: E_min - real(defReal), intent(out) :: E_max + real(defReal), intent(out) :: eMin + real(defReal), intent(out) :: eMax end subroutine energyBounds + !! + !! Make sure that trackXS of material with matIdx is at energy E = E_track + !! in ceNeutronChache + !! + !! The tracking xs correspons to the material total cross section unless TMS + !! is used. In that case, this is the material temperature majorant xs. + !! + !! Assume that call to this procedure implies that data is NOT up-to-date + !! + !! Args: + !! E [in] -> required energy [MeV] + !! matIdx [in] -> material index that needs to be updated + !! rand [inout] -> random number generator + !! + subroutine updateTrackMatXS(self, E, matIdx, rand) + import :: ceNeutronDatabase, defReal, shortInt, RNG + class(ceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + class(RNG), optional, intent(inout) :: rand + end subroutine updateTrackMatXS + !! !! Make sure that totalXS of material with matIdx is at energy E !! in ceNeutronChache @@ -196,7 +223,32 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) end subroutine updateMicroXSs !! - !! Subroutine to get the elastic scattering majorant cross section in a nuclide + !! Subroutine to retrieve the nuclide total majorant cross section over a range + !! of relative energies a nuclide can see given the material temperature. The + !! energy range is calculated based on the nuclide and material kT, on the atomic + !! mass of the nuclide, and the incident neutron energy + !! + !! The nuclide cache is updated with the 'temperature' majorant value, incident + !! neutron energy, deltakT and Doppler correction factor + !! + !! Args: + !! E [in] -> required energy [MeV] + !! kT [in] -> thermal energy of TMS material + !! nucIdx [in] -> material index that needs to be updated + !! + !! Errors: + !! FatalError if material kT is smaller than the nuclide kT + !! + subroutine updateTotalTempNucXS(self, E, kT, nucIdx) + import :: ceNeutronDatabase, defReal, shortInt + class(ceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + real(defReal), intent(in) :: kT + integer(shortInt), intent(in) :: nucIdx + end subroutine updateTotalTempNucXS + + !! + !! Function to get the elastic scattering majorant cross section in a nuclide !! over a certain energy range, defined as a function of a given temperature !! !! NOTE: This function is called by the collision operator to apply DBRC; nucIdx @@ -242,7 +294,7 @@ function getTrackingXS(self, p, matIdx, what) result(xs) select case(what) case (MATERIAL_XS) - xs = self % getTotalMatXS(p, matIdx) + xs = self % getTrackMatXS(p, matIdx) case (MAJORANT_XS) xs = self % getMajorantXS(p) @@ -252,6 +304,7 @@ function getTrackingXS(self, p, matIdx, what) result(xs) ! READ ONLY - read from previously updated cache if (p % E == trackingCache(1) % E) then xs = trackingCache(1) % xs + return else call fatalError(Here, 'Tracking cache failed to update during tracking') end if @@ -267,6 +320,37 @@ function getTrackingXS(self, p, matIdx, what) result(xs) end function getTrackingXS + !! + !! Return tracking XS for matIdx + !! + !! This is the regular material total cross section unless TMS is used. + !! If TMS is used, this is the material temperature majorant cross section. + !! + !! See nuclearDatabase_inter for details! + !! + !! Error: + !! fatalError if particle is not CE Neutron + !! + function getTrackMatXS(self, p, matIdx) result(xs) + class(ceNeutronDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + character(100),parameter :: Here = 'getTrackMatXS (ceNeutronDatabase_inter.f90)' + + ! Check dynamic type of the particle + if (p % isMG .or. p % type /= P_NEUTRON) then + call fatalError(Here, 'Dynamic type of the partcle is not CE Neutron but:'//p % typeToChar()) + end if + + ! Check Cache and update if needed + if (materialCache(matIdx) % E_track /= p % E) call self % updateTrackMatXS(p % E, matIdx, p % pRNG) + + ! Return Cross-Section + xs = materialCache(matIdx) % trackXS + + end function getTrackMatXS + !! !! Return Total XS for matIdx !! diff --git a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 index e8088fc3e..e10dc5585 100644 --- a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 +++ b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 @@ -1,9 +1,10 @@ module ceNeutronMaterial_class use numPrecision - use genericProcedures, only : fatalError - use RNG_class, only : RNG - use particle_class, only : particle + use universalVariables + use genericProcedures, only : fatalError, numToChar + use RNG_class, only : RNG + use particle_class, only : particle ! Nuclear Data Handles use materialHandle_inter, only : materialHandle @@ -12,10 +13,14 @@ module ceNeutronMaterial_class ! CE Neutron Interfaces use ceNeutronDatabase_inter, only : ceNeutronDatabase + use ceNeutronNuclide_inter, only : ceNeutronNuclide, ceNeutronNuclide_CptrCast ! Cache use ceNeutronCache_mod, only : materialCache, nuclideCache + ! Scattering procedures + use scatteringKernels_func, only : relativeEnergy_constXS, dopplerCorrectionFactor + implicit none private @@ -48,10 +53,14 @@ module ceNeutronMaterial_class !! type, public, extends(neutronMaterial) :: ceNeutronMaterial integer(shortInt) :: matIdx = 0 + real(defReal) :: kT = ZERO class(ceNeutronDatabase), pointer :: data => null() real(defReal), dimension(:), allocatable :: dens integer(shortInt), dimension(:), allocatable :: nuclides - logical(defBool) :: fissile =.false. + logical(defBool) :: fissile = .false. + logical(defBool) :: hasTMS = .false. + real(defReal) :: eUpperSab = ZERO + real(defReal) :: eLowerURR = ZERO contains ! Superclass procedures @@ -64,6 +73,7 @@ module ceNeutronMaterial_class procedure, non_overridable :: setComposition procedure, non_overridable :: getMacroXSs_byE procedure :: isFissile + procedure :: useTMS procedure, non_overridable :: sampleNuclide procedure, non_overridable :: sampleFission procedure, non_overridable :: sampleScatter @@ -80,10 +90,14 @@ elemental subroutine kill(self) class(ceNeutronMaterial), intent(inout) :: self self % matIdx = 0 + self % kT = ZERO self % data => null() - if(allocated(self % dens)) deallocate(self % dens ) - if(allocated(self % nuclides)) deallocate (self % nuclides) + if (allocated(self % dens)) deallocate(self % dens) + if (allocated(self % nuclides)) deallocate (self % nuclides) self % fissile = .false. + self % hasTMS = .false. + self % eUpperSab = ZERO + self % eLowerURR = ZERO end subroutine kill @@ -105,6 +119,7 @@ subroutine getMacroXSs_byP(self, xss, p) call fatalError(Here,'MG neutron given to CE data') end if + end subroutine getMacroXSs_byP !! @@ -118,7 +133,7 @@ end subroutine getMacroXSs_byP !! nucIdxs [in] -> correpsonding array with nucIdxs !! !! Errors: - !! FatalError if arrays have diffrent size + !! FatalError if arrays have different size !! FatalError if dens contains -ve values !! FatalError if dens has size of 0 -> no composition !! @@ -129,13 +144,13 @@ subroutine setComposition(self, dens, nucIdxs) character(100), parameter :: Here = 'setComposition (ceNeutronMaterial_class.f90)' ! Check input - if(size(dens) /= size(nucIdxs)) call fatalError(Here,'Diffrent sizes of density and nuclide vector') - if(any(dens < ZERO)) call fatalError(Here,'-ve nuclide densities are present') - if(size(dens) == 0) call fatalError(Here,'Empty composition is not allowed') + if (size(dens) /= size(nucIdxs)) call fatalError(Here,'Different sizes of density and nuclide vector') + if (any(dens < ZERO)) call fatalError(Here,'-ve nuclide densities are present') + if (size(dens) == 0) call fatalError(Here,'Empty composition is not allowed') ! Clean any current content - if(allocated(self % dens)) deallocate(self % dens) - if(allocated(self % nuclides)) deallocate(self % nuclides) + if (allocated(self % dens)) deallocate(self % dens) + if (allocated(self % nuclides)) deallocate(self % nuclides) ! Load values self % dens = dens @@ -147,25 +162,52 @@ end subroutine setComposition !! Set matIdx, pointer to a database and fissile flag !! !! All arguments are optional. Use with keyword association e.g. - !! call mat % set( matIdx = 7) + !! call mat % set(matIdx = 7) !! !! Use this procedure ONLY during build. NEVER during transport. !! IT IS NOT THREAD SAFE! !! + !! NOTE: eUpperSab and eLowerURR are fed by the aceNeutronDatabase, and they are the + !! strictest (respecitvely highest and lowest) energy limits among all nuclides + !! in the material composition. + !! !! Args: !! matIdx [in] -> material index !! database [in] -> pointer to a database that updates XSs on the ceNeutronCache - !! fissile [in] -> flag indicating whether fission data is present - !! - subroutine set(self, matIdx, database, fissile) - class(ceNeutronMaterial), intent(inout) :: self - integer(shortInt), intent(in),optional :: matIdx - class(ceNeutronDatabase),pointer, optional,intent(in) :: database - logical(defBool),intent(in), optional :: fissile + !! fissile [in] -> flag indicating whether fission data is present + !! hasTMS [in] -> flag indicating whether TMS is on + !! temp [in] -> TMS material temperature + !! eUpperSab [in] -> upper energy of S(a,b) range in the material + !! eLowerURR [in] -> lower energy of ures range in the material + !! + subroutine set(self, matIdx, database, fissile, hasTMS, temp, eUpperSab, eLowerURR) + class(ceNeutronMaterial), intent(inout) :: self + integer(shortInt), intent(in), optional :: matIdx + class(ceNeutronDatabase), pointer, optional, intent(in) :: database + logical(defBool), intent(in), optional :: fissile + logical(defBool), intent(in), optional :: hasTMS + real(defReal), intent(in), optional :: temp + real(defReal), intent(in), optional :: eUpperSab + real(defReal), intent(in), optional :: eLowerURR + character(100), parameter :: Here = 'set (ceNeutronMaterial_class.f90)' + + if (present(database)) self % data => database + if (present(fissile)) self % fissile = fissile + if (present(matIdx)) self % matIdx = matIdx + if (present(hasTMS)) self % hasTMS = hasTMS + if (present(temp)) self % kT = (kBoltzmann * temp) / joulesPerMeV + + if (present(eUpperSab)) then + if (eUpperSab < ZERO) call fatalError (Here, 'Upper Sab energy limit of material '& + &//numToChar(matIdx)//' is negative.') + self % eUpperSab = eUpperSab + end if - if(present(matIdx)) self % matIdx = matIdx - if(present(database)) self % data => database - if(present(fissile)) self % fissile = fissile + if (present(eLowerURR)) then + if (eLowerURR < ZERO) call fatalError (Here, 'Lower URR energy limit of material '& + &//numToChar(matIdx)//' is negative.') + self % eLowerURR = eLowerURR + end if end subroutine set @@ -196,13 +238,13 @@ subroutine getMacroXSs_byE(self, xss, E, rand) end subroutine getMacroXSs_byE !! - !! Return .true. if nuclide is fissile + !! Return .true. if material is fissile !! !! Args: !! None !! !! Result: - !! .TRUE. if fissile, .FALSE. otherwise + !! .true. if fissile, .false. otherwise !! !! Errors: !! None @@ -215,6 +257,28 @@ elemental function isFissile(self) result(isIt) end function isFissile + !! + !! Return .true. if TMS is on in the material and the provided energy is not + !! within the ures or S(a,b) range of any nuclide in the material. + !! + !! Args: + !! E [in] -> test energy + !! + !! Result: + !! .true. if conditions for using TMS are satisfied, .false. otherwise + !! + !! Errors: + !! None + !! + elemental function useTMS(self, E) result(shouldIt) + class(ceNeutronMaterial), intent(in) :: self + real(defReal), intent(in) :: E + logical(defBool) :: shouldIt + + shouldIt = self % hasTMS .and. (E > self % eUpperSab) .and. (E < self % eLowerURR) + + end function useTMS + !! !! Sample collision nuclide at energy E !! @@ -224,48 +288,116 @@ end function isFissile !! Args: !! E [in] -> incident energy [MeV] !! rand [inout] -> random number generator - !! - !! Result: - !! nucIdx of the sampled nuclide for collision + !! nucIdx [out] -> sampled nuclide index + !! eOut [out] -> relative energy between neutron and target (may be /= E in case of TMS) !! !! Errors: !! fatalError if sampling fails for some reason (E.G. random number > 1) !! fatalError if E is out-of-bounds of the present data !! - function sampleNuclide(self, E, rand) result(nucIdx) + subroutine sampleNuclide(self, E, rand, nucIdx, eOut) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand - integer(shortInt) :: nucIdx - real(defReal) :: xs + integer(shortInt), intent(out) :: nucIdx + real(defReal), intent(out) :: eOut + class(ceNeutronNuclide), pointer :: nuc integer(shortInt) :: i + real(defReal) :: P_acc, eMin, eMax, A, eRel, & + trackMatXS, totNucXS, dens character(100), parameter :: Here = 'sampleNuclide (ceNeutronMaterial_class.f90)' - ! Get total material XS - if(E /= materialCache(self % matIdx) % E_tot) then - call self % data % updateTotalMatXS(E, self % matIdx, rand) + ! Get material tracking XS + if (E /= materialCache(self % matIdx) % E_track) then + call self % data % updateTrackMatXS(E, self % matIdx, rand) end if - xs = materialCache(self % matIdx) % xss % total * rand % get() + trackMatXS = materialCache(self % matIdx) % trackXS * rand % get() + + ! Loop over nuclides + do i = 1,size(self % nuclides) - ! Loop over all nuclides - do i=1,size(self % nuclides) nucIdx = self % nuclides(i) - if (E /= nuclideCache(nucIdx) % E_tot) call self % data % updateTotalNucXS(E, nucIdx, rand) - xs = xs - nuclideCache(nucIdx) % xss % total * self % dens(i) - if(xs < ZERO) return + dens = self % dens(i) + + associate (nucCache => nuclideCache(nucIdx)) + + ! Retrieve nuclide XS from cache + if (self % useTMS(E)) then + + ! If the material is using TMS, the nuclide temperature majorant is needed + ! The check for the right values stored in cache happens inside the subroutine + call self % data % updateTotalTempNucXS(E, self % kT, nucIdx) + + totNucXS = nucCache % tempMajXS * nucCache % doppCorr + + else + + ! Update nuclide cache if needed + if (E /= nucCache % E_tot) call self % data % updateTotalNucXS(E, nucIdx, rand) + totNucXS = nucCache % xss % total + + end if + + trackMatXS = trackMatXS - totNucXS * dens + + ! Nuclide temporarily accepted: check TMS condition + if (trackMatXS < ZERO) then + + ! Save energy to be used to sample reaction + eOut = E + + if (self % useTMS(E)) then + + ! If the material is using TMS, retrieve nuclide and nuclide information + nuc => ceNeutronNuclide_CptrCast(self % data % getNuclide(nucIdx)) + if (.not. associated(nuc)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') + + A = nuc % getMass() + + ! Sample relative energy + eRel = relativeEnergy_constXS(E, A, nucCache % deltakT, rand) + + ! Call through system minimum and maximum energies + call self % data % energyBounds(eMin, eMax) + + ! Ensure relative energy is within energy bounds + if (eRel < eMin) eRel = eMin + if (eMax < eRel) eRel = eMax + + ! Get relative energy nuclide cross section + call self % data % updateTotalNucXS(eRel, nucIdx, rand) + + ! Calculate acceptance probability using ratio of relative energy xs to temperature majorant + P_acc = nucCache % xss % total * nucCache % doppCorr / totNucXS + + ! Accept or reject the sampled nuclide + if (rand % get() >= P_acc) nucIdx = REJECTED + + ! Overwrite energy to be used to sample reaction + eOut = eRel + + end if + + ! Exit function, return the sampled nucIdx + return + + end if + + end associate + end do ! Print error message as the inversion failed call fatalError(Here,'Nuclide sampling loop failed to terminate') - end function sampleNuclide + end subroutine sampleNuclide !! !! Sample fission nuclide given that a fission neutron was produced !! - !! Basicly samples from P(nucIdx| fission neutron produced in material) - !! Usefull when generating fission sites + !! Basically samples from P(nucIdx| fission neutron produced in material) + !! Useful when generating fission sites !! !! As such it uses nu*sigma_f !! @@ -287,30 +419,52 @@ function sampleFission(self, E, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand - integer(shortInt) :: nucIdx - real(defReal) :: xs - integer(shortInt) :: i + class(ceNeutronNuclide), pointer :: nuc + integer(shortInt) :: nucIdx, i + real(defReal) :: xs, doppCorr, A, nuckT, deltakT character(100), parameter :: Here = 'sampleFission (ceNeutronMaterial_class.f90)' ! Short-cut for nonFissile material - if(.not.self % fissile) then + if (.not. self % fissile) then nucIdx = 0 return end if ! Calculate material macroscopic nuFission - if(E /= materialCache(self % matIdx) % E_tail) then - call self % data % updateMacroXSs(E, self % matIdx, rand) - end if + ! The cache is updated without checking the energy to get the correct results with TMS + ! The relative energy flag cached is cleaned to make sure cross sections are updated + materialCache(self % matIdx) % E_rel = ZERO + call self % data % updateMacroXSs(E, self % matIdx, rand) xs = materialCache(self % matIdx) % xss % nuFission * rand % get() ! Loop over all nuclides - do i=1,size(self % nuclides) + do i = 1,size(self % nuclides) + nucIdx = self % nuclides(i) - if(E /= nuclideCache(nucIdx) % E_tail) call self % data % updateMicroXSs(E, nucIdx, rand) - xs = xs - nuclideCache(nucIdx) % xss % nuFission * self % dens(i) - if(xs < ZERO) return + + ! If the material uses TMS, the Doppler correction factor is needed + if (self % useTMS(E)) then + + nuc => ceNeutronNuclide_CptrCast(self % data % getNuclide(nucIdx)) + if (.not. associated(nuc)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') + + A = nuc % getMass() + nuckT = nuc % getkT() + deltakT = self % kT - nuckT + doppCorr = dopplerCorrectionFactor(E, A, deltakT) + + else + doppCorr = ONE + end if + + ! The nuclide cache should be at the right energy after updating the material + ! In the case of TMS where the macro xss are at a relative energy, the nuclide + ! xss to be used are at the relative energy just sampled + xs = xs - nuclideCache(nucIdx) % xss % nuFission * self % dens(i) * doppCorr + + if (xs < ZERO) return + end do ! Print error message as the inversion failed @@ -340,26 +494,49 @@ function sampleScatter(self, E, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand - integer(shortInt) :: nucIdx - real(defReal) :: xs - integer(shortInt) :: i + class(ceNeutronNuclide), pointer :: nuc + integer(shortInt) :: nucIdx, i + real(defReal) :: xs, doppCorr, A, nuckT, deltakT character(100), parameter :: Here = 'sampleScatter (ceNeutronMaterial_class.f90)' ! Calculate material macroscopic cross section of all scattering - if(E /= materialCache(self % matIdx) % E_tail) then - call self % data % updateMacroXSs(E, self % matIdx, rand) - end if + ! The cache is updated without checking the energy to get the correct results with TMS + ! The relative energy flag cached is cleaned to make sure cross sections are updated + materialCache(self % matIdx) % E_rel = ZERO + call self % data % updateMacroXSs(E, self % matIdx, rand) xs = rand % get() * (materialCache(self % matIdx) % xss % elasticScatter + & materialCache(self % matIdx) % xss % inelasticScatter) ! Loop over all nuclides - do i=1,size(self % nuclides) + do i = 1,size(self % nuclides) + nucIdx = self % nuclides(i) - if(E /= nuclideCache(nucIdx) % E_tail) call self % data % updateMicroXSs(E, nucIdx, rand) + + ! If the material uses TMS, the Doppler correction factor is needed + if (self % useTMS(E)) then + + ! If the material is using TMS, the Doppler correction factor is needed + nuc => ceNeutronNuclide_CptrCast(self % data % getNuclide(nucIdx)) + if (.not. associated(nuc)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') + + A = nuc % getMass() + nuckT = nuc % getkT() + deltakT = self % kT - nuckT + doppCorr = dopplerCorrectionFactor(E, A, deltakT) + + else + doppCorr = ONE + end if + + ! The nuclide cache should be at the right energy after updating the material + ! In the case of TMS where the macro xss are at a relative energy, the nuclide + ! xss to be used are at the relative energy just sampled xs = xs - (nuclideCache(nucIdx) % xss % elasticScatter + & - nuclideCache(nucIdx) % xss % inelasticScatter ) * self % dens(i) - if(xs < ZERO) return + nuclideCache(nucIdx) % xss % inelasticScatter ) * self % dens(i) * doppCorr + + if (xs < ZERO) return + end do ! Print error message as the inversion failed @@ -389,28 +566,51 @@ function sampleScatterWithFission(self, E, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand - integer(shortInt) :: nucIdx - real(defReal) :: xs - integer(shortInt) :: i + class(ceNeutronNuclide), pointer :: nuc + integer(shortInt) :: nucIdx, i + real(defReal) :: xs, doppCorr, A, nuckT, deltakT character(100), parameter :: Here = 'sampleScatterWithFission (ceNeutronMaterial_class.f90)' - ! Calculate material macroscopic cross section of all scattering - if(E /= materialCache(self % matIdx) % E_tail) then - call self % data % updateMacroXSs(E, self % matIdx, rand) - end if + ! Calculate material macroscopic cross section of all scattering and fission + ! The cache is updated without checking the energy to get the correct results with TMS + ! The relative energy flag cached is cleaned to make sure cross sections are updated + materialCache(self % matIdx) % E_rel = ZERO + call self % data % updateMacroXSs(E, self % matIdx, rand) xs = rand % get() * (materialCache(self % matIdx) % xss % elasticScatter + & materialCache(self % matIdx) % xss % inelasticScatter + & materialCache(self % matIdx) % xss % fission) ! Loop over all nuclides - do i=1,size(self % nuclides) + do i = 1,size(self % nuclides) + nucIdx = self % nuclides(i) - if(E /= nuclideCache(nucIdx) % E_tail) call self % data % updateMicroXSs(E, nucIdx, rand) + + ! If the material uses TMS, the Doppler correction factor is needed + if (self % useTMS(E)) then + + ! If the material is using TMS, the Doppler correction factor is needed + nuc => ceNeutronNuclide_CptrCast(self % data % getNuclide(nucIdx)) + if (.not. associated(nuc)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') + + A = nuc % getMass() + nuckT = nuc % getkT() + deltakT = self % kT - nuckT + doppCorr = dopplerCorrectionFactor(E, A, deltakT) + + else + doppCorr = ONE + end if + + ! The nuclide cache should be at the right energy after updating the material + ! In the case of TMS where the macro xss are at a relative energy, the nuclide + ! xss to be used are at the relative energy just sampled xs = xs - (nuclideCache(nucIdx) % xss % elasticScatter + & nuclideCache(nucIdx) % xss % inelasticScatter + & - nuclideCache(nucIdx) % xss % fission ) * self % dens(i) - if(xs < ZERO) return + nuclideCache(nucIdx) % xss % fission) * self % dens(i) * doppCorr + + if (xs < ZERO) return + end do ! Print error message as the inversion failed diff --git a/NuclearData/ceNeutronData/ceNeutronNuclide_inter.f90 b/NuclearData/ceNeutronData/ceNeutronNuclide_inter.f90 index 286a04e93..645c6f79e 100644 --- a/NuclearData/ceNeutronData/ceNeutronNuclide_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronNuclide_inter.f90 @@ -79,6 +79,7 @@ module ceNeutronNuclide_inter procedure(invertInelastic),deferred :: invertInelastic procedure(xsOf), deferred :: xsOf procedure(elScatteringXS), deferred :: elScatteringXS + procedure(needsSabEl), deferred :: needsSabEl end type ceNeutronNuclide @@ -153,7 +154,26 @@ function elScatteringXS(self, E) result(xs) end function elScatteringXS + !! + !! Function that checks whether this nuclide at the provided energy should + !! have S(a,b) elastic scattering data or not + !! + !! Args: + !! E [in] -> incident neutron energy + !! + !! Result: + !! True or false + !! + elemental function needsSabEl(self, E) result(doesIt) + import :: ceNeutronNuclide, defReal, defBool + class(ceNeutronNuclide), intent(in) :: self + real(defReal), intent(in) :: E + logical(defBool) :: doesIt + + end function needsSabEl + end interface + contains !! @@ -176,11 +196,11 @@ function getTotalXS(self, E, rand) result(xs) real(defReal) :: xs ! Check Cache and update if needed - if (nuclideCache(self % getNucIdx()) % E_tot /= E) then + if (nuclideCache(self % nucIdx) % E_tot /= E) then call self % data % updateTotalNucXS(E, self % nucIdx, rand) end if - xs = nuclideCache(self % getNucIdx()) % xss % total + xs = nuclideCache(self % nucIdx) % xss % total end function getTotalXS @@ -202,11 +222,11 @@ subroutine getMicroXSs(self, xss, E, rand) class(RNG), intent(inout) :: rand ! Check Cache and update if needed - if(nuclideCache(self % getNucIdx()) % E_tail /= E) then + if (nuclideCache(self % nucIdx) % E_tail /= E) then call self % data % updateMicroXSs(E, self % nucIdx, rand) end if - xss = nuclideCache(self % getNucIdx()) % xss + xss = nuclideCache(self % nucIdx) % xss end subroutine getMicroXSs @@ -240,18 +260,18 @@ subroutine set(self, nucIdx, database, fissile, mass, kT, dbrc) logical(defBool), intent(in), optional :: dbrc character(100), parameter :: Here = 'set (ceNuetronNuclide_inter.f90)' - if(present(nucIdx)) self % nucIdx = nucIdx - if(present(database)) self % data => database - if(present(fissile)) self % fissile = fissile - if(present(dbrc)) self % DBRC = dbrc + if (present(nucIdx)) self % nucIdx = nucIdx + if (present(database)) self % data => database + if (present(fissile)) self % fissile = fissile + if (present(dbrc)) self % DBRC = dbrc - if(present(mass)) then - if(mass <= ZERO) call fatalError(Here,"Mass of nuclide cannot be -ve: "//numToChar(mass)) + if (present(mass)) then + if (mass <= ZERO) call fatalError(Here,"Mass of nuclide cannot be -ve: "//numToChar(mass)) self % mass = mass end if - if(present(kT)) then - if(kT < ZERO) call fatalError(Here, "Temperature of nuclide cannot be -ve: "//numToChar(kT)) + if (present(kT)) then + if (kT < ZERO) call fatalError(Here, "Temperature of nuclide cannot be -ve: "//numToChar(kT)) self % kT = kT end if diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 729e209da..6d8460e58 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -100,6 +100,7 @@ module materialMenu_mod real(defReal),dimension(:),allocatable :: dens type(nuclideInfo),dimension(:),allocatable :: nuclides type(dictionary) :: extraInfo + logical(defBool) :: hasTMS = .false. contains procedure :: init => init_materialItem procedure :: kill => kill_materialItem @@ -200,11 +201,13 @@ subroutine display() print '(A60)', repeat('<>',30) print '(A)', "^^ MATERIAL DEFINITIONS ^^" + do i = 1,size(materialDefs) call materialDefs(i) % display() ! Print separation line print '(A)', " ><((((*> + <*))))><" end do + print '(A60)', repeat('<>',30) end subroutine display @@ -285,10 +288,20 @@ subroutine init_materialItem(self, name, idx, dict) ! Return to initial state call self % kill() - ! Load easy components c + ! Load easy components properties self % name = name self % matIdx = idx - call dict % get(self % T,'temp') + + ! Check TMS flag and read temperature + call dict % getOrDefault(self % hasTMS, 'tms', .false.) + + if (self % hasTMS .and. .not. dict % isPresent('temp')) then + call fatalError(Here, 'The material temperature must be specified when TMS is on') + end if + + call dict % getOrDefault(self % T, 'temp', ZERO) + if (self % T < ZERO) call fatalError(Here, 'The temperature of material '//numToChar(idx)//& + ' is negative: '//numToChar(self % T)) ! Get composition dictionary and load composition compDict => dict % getDictPtr('composition') @@ -298,12 +311,14 @@ subroutine init_materialItem(self, name, idx, dict) allocate(self % nuclides(size(keys))) allocate(self % dens(size(keys))) - hasSab = .false. + ! Check if S(a,b) files are specified if (dict % isPresent('moder')) then moderDict => dict % getDictPtr('moder') call moderDict % keys(moderKeys) hasSab = .true. + else + hasSab = .false. end if ! Load definitions diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index 771d0b730..082783270 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -63,6 +63,7 @@ module baseMgNeutronDatabase_class contains ! Superclass Interface procedure :: getTrackingXS + procedure :: getTrackMatXS procedure :: getTotalMatXS procedure :: getMajorantXS procedure :: matNamesMap @@ -102,7 +103,7 @@ function getTrackingXS(self, p, matIdx, what) result(xs) select case(what) case (MATERIAL_XS) - xs = self % getTotalMatXS(p, matIdx) + xs = self % getTrackMatXS(p, matIdx) case (MAJORANT_XS) xs = self % getMajorantXS(p) @@ -112,6 +113,7 @@ function getTrackingXS(self, p, matIdx, what) result(xs) ! READ ONLY - read from previously updated cache if (p % G == trackingCache(1) % G) then xs = trackingCache(1) % xs + return else call fatalError(Here, 'Tracking cache failed to update during tracking') end if @@ -128,7 +130,23 @@ function getTrackingXS(self, p, matIdx, what) result(xs) end function getTrackingXS !! - !! Get Total XS given a particle + !! Get tracking XS given a particle. In MG, it is always identical to the material + !! total XS. + !! + !! See nuclearDatabase documentation for details + !! + function getTrackMatXS(self, p, matIdx) result(xs) + class(baseMgNeutronDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + + xs = self % getTotalMatXS(p, matIdx) + + end function getTrackMatXS + + !! + !! Get total XS given a particle !! !! See nuclearDatabase documentation for details !! @@ -162,7 +180,7 @@ function getTotalMatXS(self, p, matIdx) result(xs) end function getTotalMatXS !! - !! Get Majorant XS given a particle + !! Get majorant XS given a particle !! !! See nuclearDatabase documentation for details !! diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index dea78836f..d1e25113a 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -22,6 +22,7 @@ module nuclearDatabase_inter !! !! Interface: !! getTrackingXS -> returns XS used to sample track length + !! getTrackMatXS -> returns material tracking xs, which could be different from the total (e.g., with TMS) !! getTotalMatXS -> returns total Material XS given a particle !! getMajorantXS -> returns majorant XS given particle and list of active materials !! matNamesMap -> returns pointer to map of material names to matIdx @@ -35,6 +36,7 @@ module nuclearDatabase_inter procedure(init), deferred :: init procedure(activate), deferred :: activate procedure(getTrackingXS), deferred :: getTrackingXS + procedure(getTrackMatXS), deferred :: getTrackMatXS procedure(getTotalMatXS), deferred :: getTotalMatXS procedure(getMajorantXS), deferred :: getMajorantXS procedure(matNamesMap), deferred :: matNamesMap @@ -84,7 +86,7 @@ subroutine activate(self, activeMat, silent) end subroutine activate !! - !! Return value of Tracking XS for a particle and a given request + !! Return value of tracking XS for a particle and a given request !! !! Reads all relevant state information from the particle (e.g. E or G) !! It is the XS used to sample track length: it might be the same as the @@ -112,7 +114,33 @@ function getTrackingXS(self, p, matIdx, what) result(xs) end function getTrackingXS !! - !! Return value of Material Total XS for a particle + !! Return value of materials tracking XS for a particle + !! + !! Reads all relevant state information from the particle (e.g. E or G) + !! It is the XS used to sample track length in a material: it might be the same + !! as the material total XS, or a material temperature majorant when TMS is used + !! + !! Args: + !! p [in] -> Particle at a given state + !! matIdx [in] -> Material index + !! + !! Result: + !! Value of material tracking XS [1/cm] + !! + !! Errors: + !! Undefined behaviour if the state of the particle is invalid e.g. -ve energy + !! Undefined behavior if matIdx does not correspond to a defined material + !! + function getTrackMatXS(self, p, matIdx) result(xs) + import :: nuclearDatabase, particle, shortInt, defReal + class(nuclearDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + end function getTrackMatXS + + !! + !! Return value of material total XS for a particle !! !! Reads all relevalnt state information from the particle (e.g. E or G) !! @@ -136,7 +164,7 @@ function getTotalMatXS(self, p, matIdx) result(xs) end function getTotalMatXS !! - !! Return value of Majorant XS for a particle + !! Return value of majorant XS for a particle !! !! Reads all relevalnt state information from the particle (e.g. E or G) !! Majorant XS is the largest of TRANSPORT XSs for ACTIVE materials diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index 3112f9538..e6032e73b 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -45,6 +45,7 @@ module testNeutronDatabase_class procedure :: init procedure :: activate procedure :: getTrackingXS + procedure :: getTrackMatXS procedure :: getTotalMatXS procedure :: getMajorantXS procedure :: matNamesMap @@ -179,7 +180,22 @@ function getTrackingXS(self, p, matIdx, what) result(xs) end function getTrackingXS !! - !! Return value of Material Total XS for a particle + !! Return value of material tracking XS for a particle + !! + !! See nuclearDatabase_inter for details + !! + function getTrackMatXS(self, p, matIdx) result(xs) + class(testNeutronDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + + xs = self % xsVal + + end function getTrackMatXS + + !! + !! Return value of material total XS for a particle !! !! See nuclearDatabase_inter for details !! diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index e5bd009d6..cf0f92787 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -70,6 +70,7 @@ module particle_class ! Private procedures procedure,private :: equal_particleState + end type particleState !! @@ -658,7 +659,7 @@ subroutine particleState_fromParticle(LHS,RHS) LHS % uniqueID = RHS % coords % uniqueId LHS % cellIdx = RHS % coords % cell() LHS % collisionN = RHS % collisionN - LHS % broodID = RHS % broodID + LHS % broodID = RHS % broodID end subroutine particleState_fromParticle @@ -682,7 +683,7 @@ function equal_particleState(LHS,RHS) result(isEqual) isEqual = isEqual .and. LHS % cellIdx == RHS % cellIdx isEqual = isEqual .and. LHS % uniqueID == RHS % uniqueID isEqual = isEqual .and. LHS % collisionN == RHS % collisionN - isEqual = isEqual .and. LHS % broodID == RHS % broodID + isEqual = isEqual .and. LHS % broodID == RHS % broodID if( LHS % isMG ) then isEqual = isEqual .and. LHS % G == RHS % G @@ -725,7 +726,7 @@ elemental subroutine kill_particleState(self) self % cellIdx = -1 self % uniqueID = -1 self % collisionN = 0 - self % broodID = 0 + self % broodID = 0 end subroutine kill_particleState diff --git a/SharedModules/numPrecision.f90 b/SharedModules/numPrecision.f90 index ec7d22a43..dfed315e5 100644 --- a/SharedModules/numPrecision.f90 +++ b/SharedModules/numPrecision.f90 @@ -19,13 +19,14 @@ module numPrecision ZERO = 0._defReal, & ONE = 1.0_defReal, & TWO = 2.0_defReal, & - TWO_PI = TWO * PI, & - HALF = 0.5_defReal + TWO_PI = TWO * PI, & + SQRT_PI = sqrt(PI), & + HALF = 0.5_defReal real(defReal), public, parameter :: floatTol = 1.0e-12 !*** Should be replaced real(defReal), public, parameter :: FP_REL_TOL = 1.0e-7_defReal contains - + end module numPrecision diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index b6ee050d5..1b447271c 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -64,10 +64,11 @@ module universalVariables P_NEUTRON_MG = 2 ! Search error codes - integer(shortInt), parameter :: valueOutsideArray = -1,& - tooManyIter = -2,& + integer(shortInt), parameter :: valueOutsideArray = -1, & + tooManyIter = -2, & targetNotFound = -3, & - NOT_FOUND = -3 + NOT_FOUND = -3, & + REJECTED = -4 ! Integer indexes for type of tracking cross section requested integer(shortInt), parameter :: MATERIAL_XS = 1, & @@ -78,6 +79,7 @@ module universalVariables ! Neutron mass and speed of light in vacuum from from https://physics.nist.gov/cuu/Constants/index.html real(defReal), parameter :: neutronMass = 939.56542194_defReal, & ! Neutron mass in MeV (m*c^2) lightSpeed = 2.99792458e10_defReal, & ! Light speed in cm/s + kBoltzmann = 1.380649e-23_defReal, & ! Bolztmann constant in J/K energyPerFission = 200.0_defReal ! MeV ! Unit conversion diff --git a/Tallies/TallyClerks/Tests/collisionClerk_test.f90 b/Tallies/TallyClerks/Tests/collisionClerk_test.f90 index d688c8bfd..f2086382c 100644 --- a/Tallies/TallyClerks/Tests/collisionClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/collisionClerk_test.f90 @@ -178,8 +178,9 @@ subroutine testScoring(this) call res2Dict % store('value', 1.3_defReal) ! Configure dictionary for the clerk - call clerkDict % init(6) + call clerkDict % init(7) call clerkDict % store('type','collisionClerk') + call clerkDict % store('handleVirtual', 0) call clerkDict % store(res1Name, res1Dict) call clerkDict % store(res2Name, res2Dict) @@ -297,7 +298,6 @@ subroutine testScoringVirtual(this) ! Configure dictionary for the clerk call clerkDict % init(6) call clerkDict % store('type','collisionClerk') - call clerkDict % store('handleVirtual', 1) call clerkDict % store(res1Name, res1Dict) call clerkDict % store(res2Name, res2Dict) diff --git a/Tallies/TallyClerks/Tests/keffImplicitClerk_test.f90 b/Tallies/TallyClerks/Tests/keffImplicitClerk_test.f90 index 9831fcb96..dee2593cb 100644 --- a/Tallies/TallyClerks/Tests/keffImplicitClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/keffImplicitClerk_test.f90 @@ -129,6 +129,7 @@ subroutine test1CycleBatch(this) @assertTrue(.false.,'Result is not a keffResult') end select + end subroutine test1CycleBatch !! diff --git a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 index 85a93fdd6..211039780 100644 --- a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 @@ -230,7 +230,7 @@ subroutine testScoring_clerk2(this) @assertEqual([HALF, ZERO, HALF], chi(1,:), TOL, 'Chi' ) @assertEqual([ZERO, ZERO, 4.0_defReal], transOS(1,:), TOL, 'Transport XS O.S.' ) @assertEqual([ZERO, ZERO, 5.5_defReal], transFL(1,:), TOL, 'Transport XS F.L.' ) - @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, TWO, ZERO], P0(1,:), TOL, 'P0' ) + @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, TWO, ZERO], P0(1,:), TOL, 'P0' ) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, 1.5_defReal, ZERO], P1(1,:), TOL, 'P1' ) @assertEqual([ONE, ONE, ONE, ONE, ONE, ONE, ONE, TWO, ONE], prod(1,:), TOL, 'prod' ) diff --git a/Tallies/TallyClerks/collisionClerk_class.f90 b/Tallies/TallyClerks/collisionClerk_class.f90 index 0c37c9c11..67cec49af 100644 --- a/Tallies/TallyClerks/collisionClerk_class.f90 +++ b/Tallies/TallyClerks/collisionClerk_class.f90 @@ -62,7 +62,7 @@ module collisionClerk_class integer(shortInt) :: width = 0 ! Settings - logical(defBool) :: virtual = .false. + logical(defBool) :: handleVirtual = .true. contains ! Procedures used during build @@ -120,7 +120,7 @@ subroutine init(self, dict, name) self % width = size(responseNames) ! Handle virtual collisions - call dict % getOrDefault(self % virtual,'handleVirtual', .false.) + call dict % getOrDefault(self % handleVirtual,'handleVirtual', .true.) end subroutine init @@ -134,23 +134,23 @@ elemental subroutine kill(self) call kill_super(self) ! Kill and deallocate filter - if(allocated(self % filter)) then + if (allocated(self % filter)) then deallocate(self % filter) end if ! Kill and deallocate map - if(allocated(self % map)) then + if (allocated(self % map)) then call self % map % kill() deallocate(self % map) end if ! Kill and deallocate responses - if(allocated(self % response)) then + if (allocated(self % response)) then deallocate(self % response) end if self % width = 0 - self % virtual = .false. + self % handleVirtual = .true. end subroutine kill @@ -194,29 +194,23 @@ subroutine reportInColl(self, p, xsData, mem, virtual) logical(defBool), intent(in) :: virtual type(particleState) :: state integer(shortInt) :: binIdx, i - integer(longInt) :: adrr + integer(longInt) :: addr real(defReal) :: scoreVal, flux character(100), parameter :: Here = 'reportInColl (collisionClerk_class.f90)' ! Return if collision is virtual but virtual collision handling is off - if (self % virtual) then - ! Retrieve tracking cross section from cache - flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) - else - if (virtual) return - flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) - end if + if ((.not. self % handleVirtual) .and. virtual) return ! Get current particle state state = p ! Check if within filter - if(allocated( self % filter)) then - if(self % filter % isFail(state)) return + if (allocated(self % filter)) then + if (self % filter % isFail(state)) return end if ! Find bin index - if(allocated(self % map)) then + if (allocated(self % map)) then binIdx = self % map % map(state) else binIdx = 1 @@ -225,13 +219,20 @@ subroutine reportInColl(self, p, xsData, mem, virtual) ! Return if invalid bin index if (binIdx == 0) return + ! Calculate flux with the right cross section according to virtual collision handling + if (self % handleVirtual) then + flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) + else + flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) + end if + ! Calculate bin address - adrr = self % getMemAddress() + self % width * (binIdx -1) - 1 + addr = self % getMemAddress() + self % width * (binIdx - 1) - 1 ! Append all bins - do i = 1,self % width + do i = 1, self % width scoreVal = self % response(i) % get(p, xsData) * flux - call mem % score(scoreVal, adrr + i) + call mem % score(scoreVal, addr + i) end do @@ -268,13 +269,13 @@ subroutine print(self, outFile, mem) call outFile % startBlock(self % getName()) ! If collision clerk has map print map information - if( allocated(self % map)) then + if (allocated(self % map)) then call self % map % print(outFile) end if ! Write results. ! Get shape of result array - if(allocated(self % map)) then + if (allocated(self % map)) then resArrayShape = [size(self % response), self % map % binArrayShape()] else resArrayShape = [size(self % response)] @@ -285,7 +286,7 @@ subroutine print(self, outFile, mem) call outFile % startArray(name, resArrayShape) ! Print results to the file - do i=1,product(resArrayShape) + do i = 1, product(resArrayShape) call mem % getResult(val, std, self % getMemAddress() - 1 + i) call outFile % addResult(val,std) diff --git a/Tallies/TallyClerks/collisionProbabilityClerk_class.f90 b/Tallies/TallyClerks/collisionProbabilityClerk_class.f90 index 06da8f2d8..b6182d5cb 100644 --- a/Tallies/TallyClerks/collisionProbabilityClerk_class.f90 +++ b/Tallies/TallyClerks/collisionProbabilityClerk_class.f90 @@ -21,9 +21,6 @@ module collisionProbabilityClerk_class use tallyMap_inter, only : tallyMap use tallyMapFactory_func, only : new_tallyMap - ! Tally Response - !use macroResponse_class, only : macroResponse - implicit none private @@ -72,7 +69,6 @@ module collisionProbabilityClerk_class integer(shortInt) :: N = 0 !! Number of bins !type(macroResponse) :: resp - contains ! Procedures used during build procedure :: init @@ -92,6 +88,7 @@ module collisionProbabilityClerk_class ! Deconstructor procedure :: kill + end type collisionProbabilityClerk !! @@ -192,7 +189,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) ! Find starting index in the map ! It is important that preCollision is not changed by a collisionProcessor ! before the particle is fed to the tally, otherwise results will be meaningless - sIdx = self % map % map( p % preCollision) + sIdx = self % map % map(p % preCollision) ! Find collision index in the map state = p diff --git a/Tallies/TallyClerks/keffImplicitClerk_class.f90 b/Tallies/TallyClerks/keffImplicitClerk_class.f90 index ad5f52e21..008c50e02 100644 --- a/Tallies/TallyClerks/keffImplicitClerk_class.f90 +++ b/Tallies/TallyClerks/keffImplicitClerk_class.f90 @@ -33,10 +33,9 @@ module keffImplicitClerk_class IMP_ABS = 2 ,& ! Implicit neutron absorbtion ANA_LEAK = 3 ,& ! Analog Leakage K_EFF = 4 ! k-eff estimate - !! !! A simple implicit k-eff estimator based on collison estimator of reaction rates, - !! and an analog estimators of (N,XN) reactions and leakage + !! and on analog estimators of (N,XN) reactions and leakage !! !! Private Members: !! targetSTD -> Target Standard Deviation for convergance check @@ -57,7 +56,7 @@ module keffImplicitClerk_class private real(defReal) :: targetSTD = ZERO ! Settings - logical(defBool) :: virtual = .false. + logical(defBool) :: handleVirtual = .true. contains ! Duplicate interface of the tallyClerk ! Procedures used during build @@ -78,6 +77,7 @@ module keffImplicitClerk_class procedure :: display procedure :: print procedure :: getResult + end type keffImplicitClerk contains @@ -106,7 +106,7 @@ subroutine init(self, dict, name) end if ! Handle virtual collisions - call dict % getOrDefault(self % virtual,'handleVirtual', .false.) + call dict % getOrDefault(self % handleVirtual,'handleVirtual', .true.) end subroutine init @@ -121,7 +121,7 @@ elemental subroutine kill(self) ! Kill self self % targetSTD = ZERO - self % virtual = .false. + self % handleVirtual = .true. end subroutine kill @@ -169,17 +169,18 @@ subroutine reportInColl(self, p, xsData, mem, virtual) character(100), parameter :: Here = 'reportInColl (keffImplicitClerk_class.f90)' ! Return if collision is virtual but virtual collision handling is off - if (self % virtual) then - ! Retrieve tracking cross section from cache + if ((.not. self % handleVirtual) .and. virtual) return + + ! Ensure we're not in void (could happen when scoring virtual collisions) + if (p % matIdx() == VOID_MAT) return + + ! Calculate flux with the right cross section according to virtual collision handling + if (self % handleVirtual) then flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) else - if (virtual) return flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) end if - ! Ensure we're not in void (could happen when scoring virtual collisions) - if (p % matIdx() == VOID_MAT) return - ! Get material pointer mat => neutronMaterial_CptrCast(xsData % getMaterial(p % matIdx())) if (.not.associated(mat)) then @@ -251,12 +252,12 @@ subroutine reportHist(self, p, xsData, mem) type(scoreMemory), intent(inout) :: mem real(defReal) :: histWgt - if( p % fate == leak_FATE) then + if (p % fate == leak_FATE) then ! Obtain and score history weight histWgt = p % w ! Score analog leakage - call mem % score( histWgt, self % getMemAddress() + ANA_LEAK) + call mem % score(histWgt, self % getMemAddress() + ANA_LEAK) end if @@ -274,14 +275,14 @@ subroutine reportCycleEnd(self, end, mem) integer(longInt) :: addr real(defReal) :: nuFiss, absorb, leakage, scatterMul, k_est - if( mem % lastCycle()) then + if (mem % lastCycle()) then addr = self % getMemAddress() nuFiss = mem % getScore(addr + IMP_PROD) absorb = mem % getScore(addr + IMP_ABS) leakage = mem % getScore(addr + ANA_LEAK) scatterMul = mem % getScore(addr + SCATTER_PROD) - k_est = nuFiss / (absorb + leakage - scatterMul ) + k_est = nuFiss / (absorb + leakage - scatterMul) call mem % accumulate(k_est, addr + K_EFF) end if @@ -315,7 +316,7 @@ subroutine display(self, mem) real(defReal) :: k, STD ! Get current k-eff estimate - call mem % getResult(k, STD, self % getMemAddress() + K_EFF ) + call mem % getResult(k, STD, self % getMemAddress() + K_EFF) ! Print to console print '(A,F8.5,A,F8.5)', 'k-eff (implicit): ', k, ' +/- ', STD diff --git a/Tallies/TallyClerks/mgXsClerk_class.f90 b/Tallies/TallyClerks/mgXsClerk_class.f90 index b3199f117..80ea66204 100644 --- a/Tallies/TallyClerks/mgXsClerk_class.f90 +++ b/Tallies/TallyClerks/mgXsClerk_class.f90 @@ -91,7 +91,7 @@ module mgXsClerk_class logical(defBool) :: PN = .false. ! Settings - logical(defBool) :: virtual = .false. + logical(defBool) :: handleVirtual = .true. contains ! Procedures used during build @@ -155,6 +155,9 @@ subroutine init(self, dict, name) self % width = ARRAY_SCORE_SIZE + MATRIX_SCORE_SMALL * self % energyN end if + ! Handle virtual collisions + call dict % getOrDefault(self % handleVirtual,'handleVirtual', .true.) + end subroutine init !! @@ -178,7 +181,7 @@ elemental subroutine kill(self) self % energyN = 0 self % width = 0 self % PN = .false. - self % virtual = .false. + self % handleVirtual = .false. end subroutine kill @@ -223,18 +226,12 @@ subroutine reportInColl(self, p, xsData, mem, virtual) type(neutronMacroXSs) :: xss class(neutronMaterial), pointer :: mat real(defReal) :: nuFissXS, captXS, fissXS, scattXS, flux - integer(shortInt) :: enIdx, matIdx, binIdx + integer(shortInt) :: enIdx, locIdx, binIdx integer(longInt) :: addr character(100), parameter :: Here =' reportInColl (mgXsClerk_class.f90)' ! Return if collision is virtual but virtual collision handling is off - if (self % virtual) then - ! Retrieve tracking cross section from cache - flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) - else - if (virtual) return - flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) - end if + if ((.not. self % handleVirtual) .and. virtual) return ! Get current particle state state = p @@ -248,35 +245,54 @@ subroutine reportInColl(self, p, xsData, mem, virtual) end if ! Space if (allocated(self % spaceMap)) then - matIdx = self % spaceMap % map(state) + locIdx = self % spaceMap % map(state) else - matIdx = 1 + locIdx = 1 end if ! Return if invalid bin index - if ((enIdx == self % energyN + 1) .or. matIdx == 0) return + if ((enIdx == self % energyN + 1) .or. locIdx == 0) return ! Calculate bin address - binIdx = self % energyN * (matIdx - 1) + enIdx + binIdx = self % energyN * (locIdx - 1) + enIdx addr = self % getMemAddress() + self % width * (binIdx - 1) - 1 - ! Ensure we're not in void (could happen when scoring virtual collisions) - if (p % matIdx() == VOID_MAT) return - - ! Get material pointer - mat => neutronMaterial_CptrCast(xsData % getMaterial(p % matIdx())) - if (.not.associated(mat)) then - call fatalError(Here,'Unrecognised type of material was retrived from nuclearDatabase') + ! Calculate flux with the right cross section according to virtual collision handling + if (self % handleVirtual) then + flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) + else + flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) end if - ! Retrieve material cross sections - call mat % getMacroXSs(xss, p) + ! Check if the particle is in void. This call might happen when handling virtual collisions. + ! This is relevant in the case of homogenising materials that include void: the flux + ! in void will be different than zero, and the zero reaction rates have to be averaged + if (p % matIdx() /= VOID_MAT) then - ! Calculate reaction rates - nuFissXS = xss % nuFission * flux - captXS = xss % capture * flux - fissXS = xss % fission * flux - scattXS = (xss % elasticScatter + xss % inelasticScatter) * flux + ! Get material pointer + mat => neutronMaterial_CptrCast(xsData % getMaterial(p % matIdx())) + if (.not.associated(mat)) then + call fatalError(Here,'Unrecognised type of material was retrived from nuclearDatabase') + end if + + ! Retrieve material cross sections + call mat % getMacroXSs(xss, p) + + ! Calculate reaction rates + nuFissXS = xss % nuFission * flux + captXS = xss % capture * flux + fissXS = xss % fission * flux + scattXS = (xss % elasticScatter + xss % inelasticScatter) * flux + + else + + ! Reaction rates in void are zero + nuFissXS = ZERO + captXS = ZERO + fissXS = ZERO + scattXS = ZERO + + end if ! Add scores to counters call mem % score(flux, addr + FLUX_idx) @@ -301,7 +317,7 @@ subroutine reportOutColl(self, p, MT, muL, xsData, mem) type(scoreMemory), intent(inout) :: mem type(particleState) :: preColl, postColl real(defReal) :: score, prod, mu, mu2, mu3, mu4, mu5 - integer(shortInt) :: enIdx, matIdx, binIdx, binEnOut + integer(shortInt) :: enIdx, locIdx, binIdx, binEnOut integer(longInt) :: addr ! Get pre and post collision particle state @@ -335,16 +351,16 @@ subroutine reportOutColl(self, p, MT, muL, xsData, mem) end if ! Space if (allocated(self % spaceMap)) then - matIdx = self % spaceMap % map(preColl) + locIdx = self % spaceMap % map(preColl) else - matIdx = 1 + locIdx = 1 end if ! Return if invalid bin index - if ((enIdx == self % energyN + 1) .or. matIdx == 0) return + if ((enIdx == self % energyN + 1) .or. locIdx == 0) return ! Calculate bin address - binIdx = self % energyN * (matIdx - 1) + enIdx + binIdx = self % energyN * (locIdx - 1) + enIdx addr = self % getMemAddress() + self % width * (binIdx - 1) - 1 ! Score a scattering event from group g @@ -427,7 +443,7 @@ subroutine reportSpawn(self, MT, pOld, pNew, xsData, mem) class(particleState), intent(in) :: pNew class(nuclearDatabase), intent(inout) :: xsData type(scoreMemory), intent(inout) :: mem - integer(longInt) :: addr, binIdx, enIdx, matIdx + integer(longInt) :: addr, binIdx, enIdx, locIdx if (MT == N_FISSION) then @@ -440,16 +456,16 @@ subroutine reportSpawn(self, MT, pOld, pNew, xsData, mem) end if ! Space if (allocated(self % spaceMap)) then - matIdx = self % spaceMap % map(pNew) + locIdx = self % spaceMap % map(pNew) else - matIdx = 1 + locIdx = 1 end if ! Return if invalid bin index - if ((enIdx == self % energyN + 1) .or. matIdx == 0) return + if ((enIdx == self % energyN + 1) .or. locIdx == 0) return ! Calculate bin address - binIdx = self % energyN * (matIdx - 1) + enIdx + binIdx = self % energyN * (locIdx - 1) + enIdx addr = self % getMemAddress() + self % width * (binIdx - 1) - 1 ! Score energy group of fission neutron diff --git a/Tallies/TallyClerks/simpleFMClerk_class.f90 b/Tallies/TallyClerks/simpleFMClerk_class.f90 index 93c72fe2e..a4e026aee 100644 --- a/Tallies/TallyClerks/simpleFMClerk_class.f90 +++ b/Tallies/TallyClerks/simpleFMClerk_class.f90 @@ -65,7 +65,7 @@ module simpleFMClerk_class real(defReal),dimension(:),allocatable :: startWgt integer(shortInt) :: N = 0 !! Number of bins ! Settings - logical(defBool) :: virtual = .false. + logical(defBool) :: handleVirtual = .true. contains ! Procedures used during build @@ -87,6 +87,7 @@ module simpleFMClerk_class ! Deconstructor procedure :: kill + end type simpleFMClerk !! @@ -98,7 +99,7 @@ module simpleFMClerk_class !! type,public, extends( tallyResult) :: FMresult integer(shortInt) :: N = 0 ! Size of FM - real(defReal), dimension(:,:,:),allocatable :: FM ! FM proper + real(defReal), dimension(:,:,:),allocatable :: FM ! FM proper end type FMResult contains @@ -129,7 +130,7 @@ subroutine init(self, dict, name) call self % resp % build(macroNuFission) ! Handle virtual collisions - call dict % getOrDefault(self % virtual,'handleVirtual', .false.) + call dict % getOrDefault(self % handleVirtual,'handleVirtual', .true.) end subroutine init @@ -174,11 +175,13 @@ subroutine reportCycleStart(self, start, mem) self % startWgt = ZERO ! Loop through a population and calculate starting weight in each bin - do i=1,start % popSize() - associate( state => start % get(i) ) + do i = 1,start % popSize() + + associate (state => start % get(i)) idx = self % map % map(state) - if(idx > 0) self % startWgt(idx) = self % startWgt(idx) + state % wgt + if (idx > 0) self % startWgt(idx) = self % startWgt(idx) + state % wgt end associate + end do end subroutine reportCycleStart @@ -194,21 +197,15 @@ subroutine reportInColl(self, p, xsData, mem, virtual) class(nuclearDatabase),intent(inout) :: xsData type(scoreMemory), intent(inout) :: mem logical(defBool), intent(in) :: virtual + class(neutronMaterial), pointer :: mat type(particleState) :: state integer(shortInt) :: sIdx, cIdx integer(longInt) :: addr real(defReal) :: score, flux - class(neutronMaterial), pointer :: mat - character(100), parameter :: Here = 'reportInColl simpleFMClear_class.f90' + character(100), parameter :: Here = 'reportInColl simpleFMClerk_class.f90' ! Return if collision is virtual but virtual collision handling is off - if (self % virtual) then - ! Retrieve tracking cross section from cache - flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) - else - if (virtual) return - flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) - end if + if ((.not. self % handleVirtual) .and. virtual) return ! Ensure we're not in void (could happen when scoring virtual collisions) if (p % matIdx() == VOID_MAT) return @@ -222,6 +219,13 @@ subroutine reportInColl(self, p, xsData, mem, virtual) ! Return if material is not fissile if (.not. mat % isFissile()) return + ! Calculate flux with the right cross section according to virtual collision handling + if (self % handleVirtual) then + flux = p % w / xsData % getTrackingXS(p, p % matIdx(), TRACKING_XS) + else + flux = p % w / xsData % getTotalMatXS(p, p % matIdx()) + end if + ! Find starting index in the map sIdx = self % map % map(p % preHistory) @@ -230,7 +234,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) cIdx = self % map % map(state) ! Defend against invalid collision or starting bin - if(cIdx == 0 .or. sIdx == 0 ) return + if (cIdx == 0 .or. sIdx == 0) return ! Calculate fission neutron production score = self % resp % get(p, xsData) * flux @@ -254,18 +258,18 @@ subroutine reportCycleEnd(self, end, mem) integer(longInt) :: addrFM real(defReal) :: normFactor - if(mem % lastCycle()) then + if (mem % lastCycle()) then ! Set address to the start of Fission Matrix - ! Decrease by 1 to get correct addres on the fisrt iteration of the loop + ! Decrease by 1 to get correct address on the first iteration of the loop addrFM = self % getMemAddress() - 1 ! Normalise and accumulate estimates - do i=1,self % N + do i = 1,self % N ! Calculate normalisation factor normFactor = self % startWgt(i) - if(normFactor /= ZERO) normFactor = ONE / normFactor + if (normFactor /= ZERO) normFactor = ONE / normFactor - do j=1,self % N + do j = 1,self % N ! Normalise FM column addrFM = addrFM + 1 call mem % closeBin(normFactor, addrFM) @@ -294,7 +298,7 @@ pure subroutine getResult(self, res, mem) ! Allocate result to FMresult ! Do not deallocate if already allocated to FMresult ! Its not to nice -> clean up - if(allocated(res)) then + if (allocated(res)) then select type(res) class is (FMresult) ! Do nothing @@ -315,7 +319,7 @@ pure subroutine getResult(self, res, mem) ! Check size and reallocate space if needed ! This is horrible. Hove no time to polish. Blame me (MAK) if (allocated(res % FM)) then - if( any(shape(res % FM) /= [self % N, self % N, 2])) then + if (any(shape(res % FM) /= [self % N, self % N, 2])) then deallocate(res % FM) allocate(res % FM(self % N, self % N, 2)) end if @@ -329,7 +333,7 @@ pure subroutine getResult(self, res, mem) ! Load entries addr = self % getMemAddress() - 1 do i = 1,self % N - do j=1, self % N + do j = 1, self % N addr = addr + 1 call mem % getResult(val, STD, addr) res % FM(j, i, 1) = val @@ -338,6 +342,7 @@ pure subroutine getResult(self, res, mem) end do end select + end subroutine getResult !! @@ -379,7 +384,7 @@ subroutine print(self, outFile, mem) call outFile % startArray(name, [self % N, self % N]) - do i=1,self % N * self % N + do i = 1,self % N * self % N addr = addr + 1 call mem % getResult(val, std, addr) call outFile % addResult(val, std) @@ -401,11 +406,11 @@ elemental subroutine kill(self) ! Call superclass call kill_super(self) - if(allocated(self % map)) deallocate(self % map) - if(allocated(self % startWgt)) deallocate(self % startWgt) + if (allocated(self % map)) deallocate(self % map) + if (allocated(self % startWgt)) deallocate(self % startWgt) self % N = 0 - self % virtual = .false. + self % handleVirtual = .true. call self % resp % kill() diff --git a/Tallies/TallyResponses/macroResponse_class.f90 b/Tallies/TallyResponses/macroResponse_class.f90 index a718a99af..59012772f 100644 --- a/Tallies/TallyResponses/macroResponse_class.f90 +++ b/Tallies/TallyResponses/macroResponse_class.f90 @@ -2,6 +2,7 @@ module macroResponse_class use numPrecision use endfConstants + use universalVariables, only : VOID_MAT use genericProcedures, only : fatalError, numToChar use dictionary_class, only : dictionary use particle_class, only : particle, P_NEUTRON @@ -121,14 +122,15 @@ function get(self, p, xsData) result(val) val = ZERO - ! Return 0.0 if particle is not neutron - if(p % type /= P_NEUTRON) return + ! Return zero if particle is not neutron or if the particle is in void + if (p % type /= P_NEUTRON) return + if (p % matIdx() == VOID_MAT) return ! Get pointer to active material data mat => neutronMaterial_CptrCast(xsData % getMaterial(p % matIdx())) ! Return if material is not a neutronMaterial - if(.not.associated(mat)) return + if (.not.associated(mat)) return call mat % getMacroXSs(xss, p) val = xss % get(self % MT) diff --git a/Tallies/TallyResponses/microResponse_class.f90 b/Tallies/TallyResponses/microResponse_class.f90 index 76d3d7635..225c20fce 100644 --- a/Tallies/TallyResponses/microResponse_class.f90 +++ b/Tallies/TallyResponses/microResponse_class.f90 @@ -2,6 +2,7 @@ module microResponse_class use numPrecision use endfConstants + use universalVariables, only : VOID_MAT use genericProcedures, only : fatalError, numToChar use dictionary_class, only : dictionary use particle_class, only : particle, P_NEUTRON @@ -153,17 +154,19 @@ function get(self, p, xsData) result(val) val = ZERO - ! Return 0.0 if particle is not neutron - if(p % type /= P_NEUTRON) return + ! Return zero if particle is not neutron or if the particle is in void + if (p % type /= P_NEUTRON) return + if (p % matIdx() == VOID_MAT) return ! Get pointer to active material data mat => neutronMaterial_CptrCast(xsData % getMaterial(self % matIdx)) ! Return if material is not a neutronMaterial - if(.not.associated(mat)) return + if (.not.associated(mat)) return ! Get the macroscopic cross section for the material call mat % getMacroXSs(xss, p) + ! Normalise the macroscopic cross section with the atomic density val = xss % get(self % MT) / self % dens diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 46fbb9a9e..a19f05854 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -36,6 +36,7 @@ module transportOperatorDT_class procedure :: transit => deltaTracking ! Override procedure procedure :: init + end type transportOperatorDT contains @@ -84,7 +85,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) end if ! Obtain the local cross-section - sigmaT = self % xsData % getTotalMatXS(p, p % matIdx()) + sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) ! Roll RNG to determine if the collision is real or virtual ! Exit the loop if the collision is real, report collision if virtual diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index 038a0fc1e..e7ea82435 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -33,12 +33,14 @@ module transportOperatorHT_class !! type, public, extends(transportOperator) :: transportOperatorHT real(defReal) :: cutoff ! Cutoff threshold between ST and DT + contains procedure :: transit => tracking_selection procedure, private :: deltaTracking procedure, private :: surfaceTracking ! Override procedure procedure :: init + end type transportOperatorHT contains @@ -56,7 +58,7 @@ subroutine tracking_selection(self, p, tally, thisCycle, nextCycle) majorant_inv = ONE / self % xsData % getTrackingXS(p, p % matIdx(), MAJORANT_XS) ! Obtain the local cross-section - sigmaT = self % xsData % getTotalMatXS(p, p % matIdx()) + sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) ! Calculate ratio between local cross-section and majorant ratio = sigmaT*majorant_inv @@ -114,7 +116,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) end if ! Obtain the local cross-section - sigmaT = self % xsData % getTotalMatXS(p, p % matIdx()) + sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) ! Roll RNG to determine if the collision is real or virtual ! Exit the loop if the collision is real, report collision if virtual @@ -127,6 +129,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) end do DTLoop call tally % reportTrans(p) + end subroutine deltaTracking !! diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index ddfdf5b6d..a15af04e1 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -112,7 +112,7 @@ subroutine transport(self, p, tally, thisCycle, nextCycle) call self % transit(p, tally, thisCycle, nextCycle) ! Send history reports if particle died - if( p % isDead) then + if (p % isDead) then call tally % reportHist(p) end if diff --git a/docs/User Manual.rst b/docs/User Manual.rst index d16829b97..5c0bb881e 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -776,7 +776,7 @@ Materials definition The *materials* definition is structured as: :: materials { - { temp ; + { tms <0 or 1>; temp ; composition { } *keywords* } { temp ; @@ -784,15 +784,19 @@ The *materials* definition is structured as: :: *keywords* } } -In this case, ``materialName`` can be any name chosen by the user; ``temp`` is the -material temperature in [K]. +In this case, ``materialName`` can be any name chosen by the user; the keyword ``tms`` +(*optional*, default = 0) activates Target Motion Sampling (TMS) if set to 1; TMS uses +the material temperature defined under ``temp`` [K]. ``temp`` is *optional* unless TMS +is used. .. note:: - At the moment ``temp`` is not used in any way since SCONE has no way to treat - the temperature dependence of cross-sections. It is included for future use. - To change the temperature, a user needs to set appropriate suffix to each - individual nuclide in the composition definition. + When using TMS, the temperature specified by ``temp`` must be higher than the + temperatures of the nuclides in the material composition. +.. note:: + *IMPORTANT*: When using TMS, all the tallies based on the collision estimator have to + allow scoring virtual collisions, otherwise the results will be biased. The tallies + based on the track length estimator will be biased too. The ``composition`` dictionary must always be included, but it can be empty in multi-group simulations. In continuous energy simulations, it should include a @@ -928,8 +932,12 @@ The **tally clerks** determine which kind of estimator will be used. The options that defines the domains of integration of each tally - filter (*optional*): can filter out particles with certain properties, preventing them from scoring results - - handleVirtual (*optional*, default = 0): if set to 1, delta tracking virtual collisions - are tallied with a collisionClerk as well as physical collisions + - handleVirtual (*optional*, default = 1): if set to 1, delta tracking virtual collisions + and TMS rejected collisions are tallied with a collisionClerk as well as physical collisions + +.. note:: + If TMS is on, the collisionClerk is biased for results in the TMS materials unless virtual + collisions are scored (use ) * trackClerk @@ -940,6 +948,9 @@ The **tally clerks** determine which kind of estimator will be used. The options - filter (*optional*): can filter out particles with certain properties, preventing them from scoring results +.. note:: + If TMS is on, the trackClerk is biased for results in the TMS materials + Example: :: tally { @@ -955,14 +966,18 @@ Example: :: * keffAnalogClerk, analog k_eff estimator * keffImplicitClerk, implicit k_eff estimator - - handleVirtual (*optional*, default = 0): if set to 1, delta tracking virtual collisions - are tallied with a collisionClerk as well as physical collisions + - handleVirtual (*optional*, default = 1): if set to 1, delta tracking virtual collisions + and TMS rejected collisions are tallied with a collisionClerk as well as physical collisions + +.. note:: + If TMS is on, the keffImplicitClerk is biased for results in the TMS materials unless virtual + collisions are scored (use ) Example: :: tally { k_eff1 { type keffAnalogClerk; } - k_eff2 { type keffImplicitClerk; handleVirtual 1; } + k_eff2 { type keffImplicitClerk; handleVirtual 0; } } * centreOfMassClerk, geometrical 3D center of mass estimator @@ -1009,8 +1024,12 @@ Example: :: tally map - PN (*optional*, default = 0): 1 for true; 0 for false; flag that indicates whether to calculate scattering matrices only up to P1 (``PN 0``) or P7 (``PN 1``) - - handleVirtual (*optional*, default = 0): if set to 1, delta tracking virtual collisions - are tallied with a collisionClerk as well as physical collisions + - handleVirtual (*optional*, default = 1): if set to 1, delta tracking virtual collisions + and TMS rejected collisions are tallied with a collisionClerk as well as physical collisions + +.. note:: + If TMS is on, the mgXsClerk is biased for results in the TMS materials unless virtual + collisions are scored (use ) Example: :: @@ -1039,8 +1058,12 @@ Example: :: - map: contains a dictionary with the ``tallyMap`` definition, that defines the bins of the matrix - - handleVirtual (*optional*, default = 0): if set to 1, delta tracking virtual collisions - are tallied with a collisionClerk as well as physical collisions + - handleVirtual (*optional*, default = 1): if set to 1, delta tracking virtual collisions + and TMS rejected collisions are tallied with a collisionClerk as well as physical collisions + +.. note:: + If TMS is on, the simpleFMClerk is biased for results in the TMS materials unless virtual + collisions are scored (use ) Example: ::