diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index c3db38cd2..70474ce38 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -78,6 +78,7 @@ module aceNeutronNuclide_class !! eGrid -> Energy grid for the XSs !! mainData -> Array of XSs that are required in ceNeutronMicroXSs, that is !! (total, capture, escatter, iescatter, fission, nuFission) + !! promptNuFiss -> prompt nuFission cross section for dynamic calculations !! MTdata -> array of 'reactionMT's with data for all MT reactions in the nuclide !! only reactions 1:nMT are active, that is can be sampled during tracking !! nMT -> number of active MT reactions that produce 2nd-ary neutrons @@ -110,6 +111,7 @@ module aceNeutronNuclide_class character(nameLen) :: ZAID = '' real(defReal), dimension(:), allocatable :: eGrid real(defReal), dimension(:,:), allocatable :: mainData + real(defReal), dimension(:), allocatable :: promptNuFiss type(reactionMT), dimension(:), allocatable :: MTdata integer(shortInt) :: nMT = 0 type(intMap) :: idxMT @@ -283,7 +285,7 @@ elemental subroutine kill(self) call self % elasticScatter % kill() call self % fission % kill() - if(allocated(self % MTdata)) then + if (allocated(self % MTdata)) then do i=1,size(self % MTdata) call self % MTdata(i) % kinematics % kill() end do @@ -292,9 +294,9 @@ elemental subroutine kill(self) ! Local killing self % ZAID = '' self % nMT = 0 - if(allocated(self % MTdata)) deallocate(self % MTdata) - if(allocated(self % mainData)) deallocate(self % mainData) - if(allocated(self % eGrid)) deallocate(self % eGrid) + if (allocated(self % MTdata)) deallocate(self % MTdata) + if (allocated(self % mainData)) deallocate(self % mainData) + if (allocated(self % eGrid)) deallocate(self % eGrid) call self % idxMT % kill() end subroutine kill @@ -389,10 +391,13 @@ elemental subroutine microXSs(self, xss, idx, f) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % promptNuFission = self % promptNuFiss(idx + 1) * f + (ONE-f) * self % promptNuFiss(idx) else xss % fission = ZERO xss % nuFission = ZERO + xss % promptNuFission = ZERO end if + end associate end subroutine microXSs @@ -429,9 +434,11 @@ elemental subroutine getThXSs(self, xss, idx, f, E) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % promptNuFission = self % promptNuFiss(idx + 1) * f + (ONE-f) * self % promptNuFiss(idx) else xss % fission = ZERO xss % nuFission = ZERO + xss % promptNuFission = ZERO end if ! Read S(a,b) tables for elastic scatter: return zero if eleastic scatter is off @@ -493,9 +500,11 @@ subroutine getUrrXSs(self, xss, idx, f, E, xi) if (self % isFissile()) then xss % fission = data(FISSION_XS, 2) * f + (ONE-f) * data(FISSION_XS, 1) xss % nuFission = data(NU_FISSION, 2) * f + (ONE-f) * data(NU_FISSION, 1) + xss % promptNuFission = self % promptNuFiss(idx + 1) * f + (ONE-f) * self % promptNuFiss(idx) else xss % fission = ZERO xss % nuFission = ZERO + xss % promptNuFission = ZERO end if ! Check if flag for multiplication factor (IFF) is true, and apply it to elastic scattering, @@ -525,6 +534,7 @@ subroutine getUrrXSs(self, xss, idx, f, E, xi) if(self % isFissile()) then xss % nuFission = xss % nuFission/xss % fission * val(3) + xss % promptNuFission = xss % promptNuFission/xss % fission * val(3) xss % fission = val(3) end if @@ -571,7 +581,9 @@ subroutine init(self, ACE, nucIdx, database) Ngrid = ACE % gridSize() ! Allocate space for main XSs - if(self % isFissile()) then + if (self % isFissile()) then + allocate(self % promptNuFiss(Ngrid)) + self % promptNuFiss = ZERO N = 6 else N = 4 @@ -638,6 +650,11 @@ subroutine init(self, ACE, nucIdx, database) do i = bottom, Ngrid self % mainData(NU_FISSION,i) = self % mainData(FISSION_XS,i) * & self % fission % release(self % eGrid(i)) + + ! Prompt nu fission for dynamic calculations + self % promptNuFiss(i) = self % mainData(FISSION_XS,i) * & + self % fission % releasePrompt(self % eGrid(i)) + end do end if diff --git a/NuclearData/xsPackages/neutronXsPackages_class.f90 b/NuclearData/xsPackages/neutronXsPackages_class.f90 index 196cfa4a5..89d60b0a8 100644 --- a/NuclearData/xsPackages/neutronXsPackages_class.f90 +++ b/NuclearData/xsPackages/neutronXsPackages_class.f90 @@ -35,6 +35,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: promptNuFission = ZERO contains procedure :: clean => clean_neutronMacroXSs procedure :: add => add_neutronMacroXSs @@ -62,6 +63,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: promptNuFission = ZERO contains procedure :: invert => invert_microXSs procedure :: invertBranchless => invert_microXSs_branchless @@ -89,6 +91,7 @@ elemental subroutine clean_neutronMacroXSs(self) self % capture = ZERO self % fission = ZERO self % nuFission = ZERO + self % promptNuFission = ZERO end subroutine clean_neutronMacroXSs @@ -115,6 +118,7 @@ elemental subroutine add_neutronMacroXSs(self, micro, dens) self % capture = self % capture + dens * micro % capture self % fission = self % fission + dens * micro % fission self % nuFission = self % nuFission + dens * micro % nuFission + self % promptNuFission = self % promptNuFission + dens * micro % promptNuFission end subroutine add_neutronMacroXSs @@ -145,12 +149,21 @@ elemental function get(self, MT) result(xs) case(macroEscatter) xs = self % elasticScatter + case(macroIEscatter) + xs = self % inelasticScatter + case(macroFission) xs = self % fission case(macroNuFission) xs = self % nuFission + case(macroPromptNuFission) + xs = self % promptNuFission + + case(macroDelayedNuFission) + xs = self % nuFission - self % promptNuFission + case(macroAbsorbtion) xs = self % fission + self % capture @@ -299,18 +312,24 @@ elemental function invert_microXSs_branchless(self, r) result(MT) real(defReal), intent(in):: r integer(shortInt):: MT real(defReal):: effectiveXStot - + effectiveXStot = self % elasticScatter + self % inelasticScatter & + self % nuFission - + if (r*effectiveXStot < self % elasticScatter) then MT = N_N_ELASTIC + elseif (r*effectiveXStot < self % inelasticScatter+self % elasticScatter) then MT = N_N_INELASTIC + elseif (r*effectiveXStot < self % inelasticScatter + self % elasticScatter + self % nuFission) then MT = N_FISSION - else + + else MT = huge(0_shortInt) + end if + end function invert_microXSs_branchless + end module neutronXsPackages_class diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index 2197c3274..1f30fca16 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -96,13 +96,15 @@ module endfConstants anyCapture = -201, & anyFission = -118 - ! List of Fake MT numbers for macroscopic XSs. Stolen from Serpent - integer(shortInt),parameter :: macroTotal = -1 ,& - macroCapture = -2 ,& - macroEscatter = -3 ,& - macroIEscatter = -4 ,& - macroFission = -6 ,& - macroNuFission = -7 + ! List of Fake MT numbers for macroscopic XSs. Partially stolen from Serpent + integer(shortInt),parameter :: macroTotal = -1 ,& + macroCapture = -2 ,& + macroEscatter = -3 ,& + macroIEscatter = -4 ,& + macroFission = -6 ,& + macroNuFission = -7 ,& + macroPromptNuFission = -8 ,& + macroDelayedNuFission = -9 ! List of Macro MT numbers for macroscopic XSs. Unique to SCONE (not from Serpent) integer(shortInt), parameter :: macroAllScatter = -20 ,& diff --git a/Tallies/TallyResponses/macroResponse_class.f90 b/Tallies/TallyResponses/macroResponse_class.f90 index a718a99af..e83ddb9b1 100644 --- a/Tallies/TallyResponses/macroResponse_class.f90 +++ b/Tallies/TallyResponses/macroResponse_class.f90 @@ -88,14 +88,13 @@ subroutine build(self, MT) ! Check that MT number is valid select case(MT) - case(macroTotal, macroCapture, macroFission, macroNuFission, macroAbsorbtion) + case(macroTotal, macroCapture, macroFission, macroNuFission, macroAbsorbtion, & + macroEscatter, macroIEscatter, macroPromptNuFission, macroDelayedNuFission) ! Do nothing. MT is Valid - case(macroEscatter) - call fatalError(Here,'Macroscopic Elastic scattering is not implemented yet') - case default call fatalError(Here,'Unrecognised MT number: '// numToChar(self % MT)) + end select ! Load MT @@ -122,15 +121,16 @@ function get(self, p, xsData) result(val) val = ZERO ! Return 0.0 if particle is not neutron - if(p % type /= P_NEUTRON) return + if (p % type /= P_NEUTRON) 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) end function get