diff --git a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 index 3d51ee384..ed7bd9d9c 100644 --- a/NuclearData/Reactions/reactionMG/fissionMG_class.f90 +++ b/NuclearData/Reactions/reactionMG/fissionMG_class.f90 @@ -18,6 +18,8 @@ module fissionMG_class !! public :: fissionMG_TptrCast + real(defReal), private, parameter :: KAPPA_DEFAULT = 202.27 ! [MeV] + !! !! A special type of MG reaction that contains data related to fission !! @@ -26,14 +28,17 @@ module fissionMG_class !! Fission spectrum is independent of the energy group !! !! Public members: - !! data -> data for fissionMG [energyGroup, reaction (nu or chi)] + !! data -> data for fissionMG [energyGroup, reaction (nu or chi)] + !! kappa -> Energy per fission. 200 MeV by default. !! !! Interface: !! reactionMG interface !! buildFromDict -> builds fissionMG from a SCONE dictionary + !! getKappa -> obtains the energy per fission in MeV !! type, public, extends(reactionMG) :: fissionMG real(defReal),dimension(:,:),allocatable :: data + real(defReal) :: kappa = KAPPA_DEFAULT contains ! Superclass procedures procedure :: init @@ -46,6 +51,7 @@ module fissionMG_class ! Local procedures procedure :: buildFromDict + procedure :: getKappa end type fissionMG @@ -216,6 +222,9 @@ subroutine buildFromDict(self, dict) ! Get number of groups call dict % get(nG, 'numberOfGroups') + + ! Get energy per fission + call dict % getOrDefault(self % kappa, 'kappa', KAPPA_DEFAULT) ! Allocate space allocate(self % data(nG, 2)) @@ -246,6 +255,18 @@ subroutine buildFromDict(self, dict) end subroutine buildFromDict + !! + !! Get the energy release from fission + !! + pure function getKappa(self) result(kappa) + class(fissionMG), intent(in) :: self + real(defReal) :: kappa + + kappa = self % kappa + + end function getKappa + + !! !! Cast reactionHandle pointer to fissionMG pointer !! diff --git a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 index bb64e01a2..c782f9672 100644 --- a/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 +++ b/NuclearData/Reactions/uncorrelatedReactionCE/fissionCE_class.f90 @@ -54,6 +54,7 @@ module fissionCE_class !! eLawPrompt -> Energy Law for the prompt fission neutrons !! nuBarDelayed -> Delayed release table for incindent energy [MeV] !! delayed -> Information about Delayed emission precursors + !! Q -> Q value for the reaction [MeV] !! !! Interface: !! uncorrelatedReactionCE interface @@ -65,6 +66,7 @@ module fissionCE_class class(energyLawENDF), allocatable :: eLawPrompt class(releaseLawENDF),allocatable :: nuBarDelayed type(precursor),dimension(:),allocatable :: delayed + real(defReal) :: Q = ZERO contains ! Superclass procedures @@ -79,6 +81,7 @@ module fissionCE_class ! Type specific procedures procedure :: buildFromACE + procedure :: getQ end type fissionCE contains @@ -148,8 +151,21 @@ elemental subroutine kill(self) deallocate(self % delayed) end if + self % Q = ZERO + end subroutine kill + !! + !! Returns the Q-value + !! + pure function getQ(self) result(Q) + class(fissionCE), intent(in) :: self + real(defReal) :: Q + + Q = self % Q + + end function getQ + !! !! Returns true if reaction is in Centre-Of-Mass frame !! @@ -349,6 +365,7 @@ subroutine buildFromACE(self, ACE, MT) ! Read basic data call new_totalNU(self % nuBarTotal, ACE) call new_energyLawENDF(self % eLawPrompt, ACE, MT) + self % Q = ACE % QforMT(MT) ! Read Delayed Data if (withDelayed) then diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index e503fc6a5..8492876d5 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -937,6 +937,9 @@ subroutine init(self, dict, ptr, silent ) end if + ! Obtain the energy per fission with which to scale fission heating + call dict % getOrDefault(self % H235, 'energyPerFission', 202.27_defReal) + ! Get path to ACE library call dict % get(aceLibPath,'aceLibrary') diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index d94e428d3..6486c987c 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -41,14 +41,17 @@ module aceNeutronNuclide_class ! Grid location parameters - integer(shortInt), parameter :: TOTAL_XS = 1 - integer(shortInt), parameter :: ESCATTER_XS = 2 - integer(shortInt), parameter :: IESCATTER_XS = 3 - integer(shortInt), parameter :: CAPTURE_XS = 4 - integer(shortInt), parameter :: FISSION_XS = 5 - integer(shortInt), parameter :: NU_FISSION = 6 - integer(shortInt), parameter :: PROMPT_NU_FISSION = 7 + integer(shortInt), parameter :: TOTAL_XS = 1 + integer(shortInt), parameter :: ESCATTER_XS = 2 + integer(shortInt), parameter :: IESCATTER_XS = 3 + integer(shortInt), parameter :: CAPTURE_XS = 4 + integer(shortInt), parameter :: FISSION_XS = 5 + integer(shortInt), parameter :: NU_FISSION = 6 + integer(shortInt), parameter :: KAPPA_XS = 7 + integer(shortInt), parameter :: PROMPT_NU_FISSION = 8 + integer(shortInt), parameter :: NON_FISSILE_SIZE = 4, & + FISSILE_SIZE = 8 !! !! Groups data related to an MT reaction @@ -437,10 +440,12 @@ 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 % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) xss % promptNuFission = data(PROMPT_NU_FISSION, 2) * f + (ONE-f) * data(PROMPT_NU_FISSION, 1) else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO xss % promptNuFission = ZERO end if end associate @@ -486,10 +491,12 @@ subroutine getThXSs(self, xss, idx, f, E, kT, rand) 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 % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) xss % promptNuFission = data(PROMPT_NU_FISSION, 2) * f + (ONE-f) * data(PROMPT_NU_FISSION, 1) else - xss % fission = ZERO - xss % nuFission = ZERO + xss % fission = ZERO + xss % nuFission = ZERO + xss % kappaXS = ZERO xss % promptNuFission = ZERO end if @@ -560,10 +567,12 @@ 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 % kappaXS = data(KAPPA_XS, 2) * f + (ONE-f) * data(KAPPA_XS, 1) xss % promptNuFission = data(PROMPT_NU_FISSION, 2) * f + (ONE-f) * data(PROMPT_NU_FISSION, 1) else - xss % fission = ZERO - xss % nuFission = ZERO + xss % fission = ZERO + xss % nuFission = ZERO + xss % kappaXS = ZERO xss % promptNuFission = ZERO end if @@ -594,6 +603,7 @@ subroutine getUrrXSs(self, xss, idx, f, E, xi) if (self % isFissile()) then xss % nuFission = xss % nuFission/xss % fission * val(3) + xss % kappaXS = xss % kappaXS / xss % fission * val(3) xss % promptNuFission = xss % promptNuFission/xss % fission * val(3) xss % fission = val(3) end if @@ -733,6 +743,7 @@ subroutine init(self, ACE, nucIdx, database) top, firstIdxMT4 real(defReal), dimension(:), allocatable :: xsMT4 type(stackInt) :: scatterMT, absMT + real(defReal) :: H_Q character(100), parameter :: Here = "init (aceNeutronNuclide_class.f90)" ! Reset nuclide just in case @@ -753,9 +764,9 @@ subroutine init(self, ACE, nucIdx, database) ! Allocate space for main XSs if (self % isFissile()) then - N = 7 + N = FISSILE_SIZE else - N = 4 + N = NON_FISSILE_SIZE end if allocate(self % mainData(N, Ngrid)) @@ -815,10 +826,20 @@ subroutine init(self, ACE, nucIdx, database) call self % fission % init(ACE, N_f) end if - ! Calculate nuFission + ! Obtain Heating/Q scaling ratio + ! Check if database is associated in order to satisfy tests where it might not be! + if (associated(database)) then + H_Q = database % H235 / database % Q235 + else + H_Q = ONE + end if + + ! Calculate nuFission and kappaXS do i = bottom, Ngrid - self % mainData(NU_FISSION,i) = self % mainData(FISSION_XS,i) * & - self % fission % release(self % eGrid(i)) + self % mainData(NU_FISSION,i) = self % mainData(FISSION_XS,i) * & + self % fission % release(self % eGrid(i)) + self % mainData(KAPPA_XS,i) = self % mainData(FISSION_XS,i) * & + self % fission % getQ() * H_Q self % mainData(PROMPT_NU_FISSION,i) = self % mainData(FISSION_XS,i) * & self % fission % releasePrompt(self % eGrid(i)) end do diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index d9352adb9..96d2c41fd 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -39,6 +39,8 @@ module ceNeutronDatabase_inter !! !! Public Members: !! mapDBRCnuc -> map to link indexes of DBRC nuclides with their corresponding 0K + !! H235 -> heating fission of U235, used for scaling fission energy + !! Q235 -> Q-value of U235 fission, used for scaling fission energy !! !! Interface: !! nuclearDatabase Interface @@ -53,6 +55,10 @@ module ceNeutronDatabase_inter !! type, public, abstract, extends(nuclearDatabase) :: ceNeutronDatabase type(intMap) :: mapDBRCnuc + + ! Default fission scaling [MeV] + real(defReal) :: H235 = 202.27_defReal + real(defReal) :: Q235 = 193.406_defReal contains diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 index dc98bdc35..068a16549 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronMaterial_class.f90 @@ -129,9 +129,11 @@ subroutine getMacroXSs_byG(self, xss, G, rand) if(self % isFissile()) then xss % fission = self % data(FISSION_XS, G) xss % nuFission = self % data(NU_FISSION, G) + xss % kappaXS = xss % fission * self % fission % getKappa() else xss % fission = ZERO xss % nuFission = ZERO + xss % kappaXS = ZERO end if end subroutine getMacroXSs_byG diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index 32f17eb62..ee448b9f7 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -71,11 +71,12 @@ module testNeutronDatabase_class !! captureXS [in] -> Optional. Value of Capture XS !! fissionXS [in] -> Optional. Value of Fission XS !! nuFissionXS [in] -> Optional Value of nuFission + !! kappaXS [in] -> Optional Value of kappa XS !! !! Errors: !! None !! - subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuFissionXS) + subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuFissionXS, kappaXS) class(testNeutroNDatabase), intent(inout) :: self real(defReal), intent(in) :: xsVal real(defReal), intent(in),optional :: eScatterXS @@ -83,6 +84,7 @@ subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuF real(defReal), intent(in),optional :: captureXS real(defReal), intent(in),optional :: fissionXS real(defReal), intent(in),optional :: nuFissionXS + real(defReal), intent(in),optional :: kappaXS self % xsVal = xsVal @@ -126,6 +128,13 @@ subroutine build(self, xsVal, eScatterXS, ieScatterXS ,captureXS, fissionXS, nuF else self % mat % xss % nuFission = xsVal end if + + ! kappa * fission + if(present(kappaXS)) then + self % mat % xss % kappaXS = kappaXS + else + self % mat % xss % kappaXS = xsVal + end if end subroutine build diff --git a/NuclearData/xsPackages/neutronXsPackages_class.f90 b/NuclearData/xsPackages/neutronXsPackages_class.f90 index 1ba424c6c..f48ee9d0e 100644 --- a/NuclearData/xsPackages/neutronXsPackages_class.f90 +++ b/NuclearData/xsPackages/neutronXsPackages_class.f90 @@ -2,7 +2,6 @@ !! This module brakes standard rules !! It contains a library of XS Packages for Neutron particle type !! -!! module neutronXsPackages_class use numPrecision @@ -22,6 +21,7 @@ module neutronXsPackages_class !! capture -> sum of all reactions without secendary neutrons excluding fission [1/cm] !! fission -> total Fission MT=18 Cross-section [1/cm] !! nuFission -> total average neutron production Cross-section [1/cm] + !! kappaXS -> heating cross-section [MeV/cm] !! promptNuFission -> prompt-only average neutron production Cross-section [1/cm] !! !! Interface: @@ -36,6 +36,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: kappaXS = ZERO real(defReal) :: promptNuFission = ZERO contains procedure :: clean => clean_neutronMacroXSs @@ -56,6 +57,7 @@ module neutronXsPackages_class !! capture -> all reactions without secendary neutrons excluding fission [barn] !! fission -> total Fission MT=18 Cross-section [barn] !! nuFission -> total average neutron production Cross-section [barn] + !! kappaXS -> heating Cross-section [MeV*barn] !! promptNuFission -> prompt-only average neutron production Cross-section [barn] !! type, public :: neutronMicroXSs @@ -65,6 +67,7 @@ module neutronXsPackages_class real(defReal) :: capture = ZERO real(defReal) :: fission = ZERO real(defReal) :: nuFission = ZERO + real(defReal) :: kappaXS = ZERO real(defReal) :: promptNuFission = ZERO contains procedure :: invert => invert_microXSs @@ -92,6 +95,7 @@ elemental subroutine clean_neutronMacroXSs(self) self % capture = ZERO self % fission = ZERO self % nuFission = ZERO + self % kappaXS = ZERO self % promptNuFission = ZERO end subroutine clean_neutronMacroXSs @@ -119,6 +123,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 % kappaXS = self % kappaXS + dens * micro % kappaXS self % promptNuFission = self % promptNuFission + dens * micro % promptNuFission end subroutine add_neutronMacroXSs @@ -165,6 +170,9 @@ elemental function get(self, MT) result(xs) case(macroNuFission) xs = self % nuFission + case(macroKappaFission) + xs = self % kappaXS + case(macroPromptNuFission) xs = self % promptNuFission @@ -245,7 +253,7 @@ end function invert_macroXSs !! !! Use a real r in <0;1> to sample reaction from Microscopic XSs !! - !! This function involves a bit of code so is written for conviniance + !! This function involves a bit of code so is written for convenience !! !! Args: !! r [in] -> Real number in <0;1> diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index e65217102..de8ee3dba 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -87,6 +87,7 @@ module endfConstants N_Xt = 205 ,& N_X3He = 206 ,& N_Xa = 207 ,& + N_KAPPA = 301 ,& nubar_tot = 452 ,& nubar_del = 455 ,& nubar_prompt = 456 ,& @@ -119,7 +120,8 @@ module endfConstants macroFission = -6 ,& macroNuFission = -7 ,& macroPromptNuFission = -8 ,& - macroDelayedNuFission = -9 + macroDelayedNuFission = -9 ,& + macroKappaFission = -80 ! List of Macro MT numbers for macroscopic XSs. Unique to SCONE (not from Serpent) integer(shortInt), parameter :: macroAllScatter = -20 ,& @@ -134,10 +136,10 @@ module endfConstants 36, 37, 38, 41, 42, 44, 45, 91, 101, 102, 103, & 104, 105, 106, 107, 108, 109, 111, 112, 113, & 114, 115, 116, 117, 203, 204, 205, 206, 207, & - 891, N_Nl, N_2Nl] + 301, 891, N_Nl, N_2Nl] integer(shortInt), dimension(*), parameter :: availableMacroMTs = & - [-1, -2, -3, -4, -6, -7, -20, -21, -22] + [-1, -2, -3, -4, -6, -7, -20, -21, -22, -80] diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 7986d4fc3..5a60419ab 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -85,8 +85,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 + kBoltzmann = 1.380649e-23_defReal, & ! Boltzmann constant in J/K kBoltzmannMeV = kBoltzmann / joulesPerMeV ! Global name variables used to define specific geometry or field types diff --git a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 index d14e284c4..d193cdb13 100644 --- a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 @@ -74,7 +74,7 @@ subroutine setUp(this) ! Build test neutronDatabase call this % nucData % build(ONE, captureXS = 2.0_defReal, & - fissionXS = 1.5_defReal, nuFissionXS = 3.0_defReal) + fissionXS = 1.5_defReal, nuFissionXS = 3.0_defReal, kappaXS = 300.0_defReal) end subroutine setUp @@ -105,8 +105,8 @@ subroutine testScoring_clerk1(this) type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, P2, P3, P4, & - P5, P6, P7, prod + nu, chi, kappa, P0, P1, P2, P3, & + P4, P5, P6, P7, prod real(defReal), parameter :: TOL = 1.0E-9 ! Configure memory @@ -143,7 +143,7 @@ subroutine testScoring_clerk1(this) call mem % closeCycle(ONE) ! Process and get results - call this % clerk_test1 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call this % clerk_test1 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) call this % clerk_test1 % processPN(mem, P2, P3, P4, P5, P6, P7) ! Verify results of scoring @@ -151,6 +151,7 @@ subroutine testScoring_clerk1(this) @assertEqual([ZERO, ZERO, 1.5_defReal, ZERO], fiss(1,:), TOL, 'Fission XS' ) @assertEqual([ZERO, ZERO, TWO, ZERO], nu(1,:), TOL, 'NuFission XS' ) @assertEqual([ZERO, ZERO, HALF, HALF], chi(1,:), TOL, 'Chi' ) + @assertEqual([ZERO, ZERO, 300.0_defReal, ZERO], kappa(1,:), TOL, 'Kappa XS' ) @assertEqual([ZERO, ZERO, 4.0_defReal, ZERO], transOS(1,:), TOL, 'Transport XS O.S.' ) @assertEqual([ZERO, ZERO, 5.5_defReal, ZERO], transFL(1,:), TOL, 'Transport XS F.L.' ) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, TWO, ZERO, ZERO], P0(1,:), TOL, 'P0' ) @@ -165,7 +166,7 @@ subroutine testScoring_clerk1(this) @assertEqual([ZERO, ZERO, ZERO, ZERO, ZERO, -0.0683670044_defReal, ZERO, ZERO], P7(1,:), TOL, 'P7' ) ! Test getting size - @assertEqual(100, this % clerk_test1 % getSize(),'Test getSize():') + @assertEqual(104, this % clerk_test1 % getSize(),'Test getSize():') ! Test correctness of output calls call out % init('dummyPrinter', fatalErrors = .false.) @@ -185,7 +186,7 @@ subroutine testScoring_clerk2(this) type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, prod + nu, chi, kappa, P0, P1, prod real(defReal), parameter :: TOL = 1.0E-9 ! Configure memory @@ -221,13 +222,14 @@ subroutine testScoring_clerk2(this) call mem % closeCycle(ONE) ! Process and get results - call this % clerk_test2 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call this % clerk_test2 % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! Verify results of scoring @assertEqual([ZERO, ZERO, TWO], capt(1,:), TOL, 'Capture XS' ) @assertEqual([ZERO, ZERO, 1.5_defReal], fiss(1,:), TOL, 'Fission XS' ) @assertEqual([ZERO, ZERO, TWO], nu(1,:), TOL, 'NuFission XS' ) @assertEqual([HALF, ZERO, HALF], chi(1,:), TOL, 'Chi' ) + @assertEqual([ZERO, ZERO, 300.0_defReal], kappa(1,:), TOL, 'Kappa XS' ) @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' ) @@ -235,7 +237,7 @@ subroutine testScoring_clerk2(this) @assertEqual([ONE, ONE, ONE, ONE, ONE, ONE, ONE, TWO, ONE], prod(1,:), TOL, 'prod' ) ! Test getting size - @assertEqual(48, this % clerk_test2 % getSize(),'Test getSize():') + @assertEqual(51, this % clerk_test2 % getSize(),'Test getSize():') ! Test correctness of output calls call out % init('dummyPrinter', fatalErrors = .false.) diff --git a/Tallies/TallyClerks/mgXsClerk_class.f90 b/Tallies/TallyClerks/mgXsClerk_class.f90 index a65d4aaf0..567cbc1e6 100644 --- a/Tallies/TallyClerks/mgXsClerk_class.f90 +++ b/Tallies/TallyClerks/mgXsClerk_class.f90 @@ -29,7 +29,7 @@ module mgXsClerk_class private !! Size of clerk memory - integer(shortInt), parameter :: ARRAY_SCORE_SIZE = 7 ,& ! Size of data to store as 1D arrays + integer(shortInt), parameter :: ARRAY_SCORE_SIZE = 8 ,& ! Size of data to store as 1D arrays MATRIX_SCORE_SMALL = 3 ,& ! Size of data to store as 2D arrays when scoring up to P1 MATRIX_SCORE_FULL = 9 ! Size of data to store as 2D arrays when scoring up to P7 @@ -40,14 +40,16 @@ module mgXsClerk_class FISS_idx = 4 ,& ! Fission macroscopic reaction rate NUBAR_idx = 5 ,& ! NuBar CHI_idx = 6 ,& ! Fission neutron spectrum - SCATT_EV_idx = 7 ! Analog: number of scattering events + KAPPAFISS_idx = 7 ,& ! Fission heating macroscopic reaction rate + SCATT_EV_idx = 8 ! Analog: number of scattering events !! !! Multi-group macroscopic cross section calculation !! !! It prints out: - !! capture xs, fission xs, transport xs, nu, chi, the P0 and P1 scattering matrices, - !! the P0 scattering production matrix. On request, also the P2 -> P7 scattering matrices. + !! capture xs, fission xs, transport xs, nu, chi, kappa * fission xs, the P0 and P1 scattering + !! matrices, the P0 scattering production matrix. On request, also the P2 -> P7 scattering + !! matrices. !! !! NOTE: !! - the cross sections are tallied with a collision estimator; @@ -230,7 +232,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) type(particleState) :: state type(neutronMacroXSs) :: xss class(neutronMaterial), pointer :: mat - real(defReal) :: nuFissXS, captXS, fissXS, scattXS, flux + real(defReal) :: nuFissXS, captXS, fissXS, scattXS, kappaXS, flux integer(shortInt) :: enIdx, locIdx, binIdx integer(longInt) :: addr character(100), parameter :: Here =' reportInColl (mgXsClerk_class.f90)' @@ -288,6 +290,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) captXS = xss % capture * flux fissXS = xss % fission * flux scattXS = (xss % elasticScatter + xss % inelasticScatter) * flux + kappaXS = xss % kappaXS * flux else @@ -296,6 +299,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) captXS = ZERO fissXS = ZERO scattXS = ZERO + kappaXS = ZERO end if @@ -305,6 +309,7 @@ subroutine reportInColl(self, p, xsData, mem, virtual) call mem % score(captXS, addr + CAPT_idx) call mem % score(fissXS, addr + FISS_idx) call mem % score(scattXS, addr + SCATT_idx) + call mem % score(kappaXS, addr + KAPPAFISS_idx) end subroutine reportInColl @@ -500,7 +505,7 @@ end subroutine reportSpawn !! none !! pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_res, & - nu_res, chi_res, P0_res, P1_res, prod_res) + nu_res, chi_res, kappa_res, P0_res, P1_res, prod_res) class(mgXsClerk), intent(in) :: self type(scoreMemory), intent(in) :: mem real(defReal), dimension(:,:), allocatable, intent(out) :: capt_res @@ -509,6 +514,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r real(defReal), dimension(:,:), allocatable, intent(out) :: transOS_res real(defReal), dimension(:,:), allocatable, intent(out) :: nu_res real(defReal), dimension(:,:), allocatable, intent(out) :: chi_res + real(defReal), dimension(:,:), allocatable, intent(out) :: kappa_res real(defReal), dimension(:,:), allocatable, intent(out) :: P0_res real(defReal), dimension(:,:), allocatable, intent(out) :: P1_res real(defReal), dimension(:,:), allocatable, intent(out) :: prod_res @@ -518,7 +524,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r integer(shortInt) :: N, M, i, j, k, g1, gEnd, idx real(defReal) :: capt, fiss, scatt, nu, chi, P0, P1, prod, sumChi, flux, scattProb, & captStd, fissStd, scattStd, nuStd, chiStd, P0std, P1std, prodStd, & - fluxStd, scattProbStd, scattXS, scattXSstd + fluxStd, scattProbStd, scattXS, scattXSstd, kappa, kappaStd ! Get number of bins N = self % energyN @@ -526,9 +532,9 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r ! Allocate arrays for MG xss allocate( capt_res(2,N*M), fiss_res(2,N*M), transFL_res(2,N*M), transOS_res(2,N*M), & - nu_res(2,N*M), chi_res(2,N*M), P0_res(2,N*N*M), P1_res(2,N*N*M), & - prod_res(2,N*N*M), tot(N*M), fluxG(N*M), delta(M, N), totStd(N*M), & - fluxGstd(N*M), deltaStd(M, N) ) + nu_res(2,N*M), chi_res(2,N*M), kappa_res(2,N*M), P0_res(2,N*N*M), & + P1_res(2,N*N*M), prod_res(2,N*N*M), tot(N*M), fluxG(N*M), delta(M, N), & + totStd(N*M), fluxGstd(N*M), deltaStd(M, N) ) ! Initialise values sumChi = 0 ! to normalise chi @@ -548,6 +554,7 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r call mem % getResult(nu, nuStd, addr + NUBAR_idx) call mem % getResult(chi, chiStd, addr + CHI_idx) call mem % getResult(scattProb, scattProbStd, addr + SCATT_EV_idx) + call mem % getResult(kappa, kappaStd, addr + KAPPAFISS_idx) ! Calculate MG constants, being careful to avoid division by zero ! If flux is zero all cross sections must be set to zero @@ -564,13 +571,16 @@ pure subroutine processRes(self, mem, capt_res, fiss_res, transFL_res, transOS_r ! Calculate fission production term and uncertainties if (fiss == ZERO) then - fiss_res(1:2,i) = ZERO - nu_res(1:2,i) = ZERO + fiss_res(1:2,i) = ZERO + nu_res(1:2,i) = ZERO + kappa_res(1:2,i) = ZERO else fiss_res(1,i) = fiss/flux fiss_res(2,i) = fiss_res(1,i) * sqrt((fissStd/fiss)**2 + (fluxStd/flux)**2) nu_res(1,i) = nu/fiss nu_res(2,i) = nu_res(1,i) * sqrt((nuStd/nu)**2 + (fissStd/fiss)**2) + kappa_res(1,i) = kappa/flux + kappa_res(2,i) = kappa_res(1,i) * sqrt((kappaStd/kappa)**2 + (fluxStd/flux)**2) end if ! Store fission spectrum @@ -758,8 +768,8 @@ subroutine print(self, outFile, mem) type(scoreMemory), intent(in) :: mem integer(shortInt),dimension(:),allocatable :: resArrayShape real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & - nu, chi, P0, P1, P2, P3, P4, & - P5, P6, P7, prod + nu, chi, kappa, P0, P1, P2, P3, & + P4, P5, P6, P7, prod character(nameLen) :: name integer(shortInt) :: i @@ -784,7 +794,7 @@ subroutine print(self, outFile, mem) end if ! Process and get results - call self % processRes(mem, capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + call self % processRes(mem, capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! Print results name = 'capture' @@ -801,6 +811,13 @@ subroutine print(self, outFile, mem) end do call outFile % endArray() + name = 'kappaFission' + call outFile % startArray(name, resArrayShape) + do i=1,product(resArrayShape) + call outFile % addResult(kappa(1,i),kappa(2,i)) + end do + call outFile % endArray() + name = 'transportFluxLimited' call outFile % startArray(name, resArrayShape) do i=1,product(resArrayShape) @@ -853,7 +870,7 @@ subroutine print(self, outFile, mem) call outFile % endArray() ! Deallocate to limit memory consumption when writing to the output file - deallocate(capt, fiss, transFL, transOS, nu, chi, P0, P1, prod) + deallocate(capt, fiss, transFL, transOS, nu, chi, kappa, P0, P1, prod) ! If high order scattering is requested, print the other matrices if (self % PN) then diff --git a/Tallies/TallyResponses/Tests/macroResponse_test.f90 b/Tallies/TallyResponses/Tests/macroResponse_test.f90 index bfa87292e..69c3815f0 100644 --- a/Tallies/TallyResponses/Tests/macroResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/macroResponse_test.f90 @@ -18,6 +18,7 @@ module macroResponse_test type(macroResponse) :: response_fission type(macroResponse) :: response_nuFission type(macroResponse) :: response_absorbtion + type(macroResponse) :: response_kappa type(macroResponse) :: response_MT type(testNeutronDatabase) :: xsData contains @@ -36,8 +37,8 @@ subroutine setUp(this) type(dictionary) :: tempDict ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission - call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal) + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappaFission + call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set up responses ! Total @@ -75,6 +76,13 @@ subroutine setUp(this) call this % response_absorbtion % init(tempDict) call tempDict % kill() + ! Energy deposition + call tempDict % init(2) + call tempDict % store('type','macroResponse') + call tempDict % store('MT', macroKappaFission) + call this % response_kappa % init(tempDict) + call tempDict % kill() + ! MT number call tempDict % init(2) call tempDict % store('type','macroResponse') @@ -117,6 +125,7 @@ subroutine testGettingResponse(this) @assertEqual(1.0_defReal, this % response_fission % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_nuFission % get(p, this % xsData), TOL) @assertEqual(3.0_defReal, this % response_absorbtion % get(p, this % xsData), TOL) + @assertEqual(9.0_defReal, this % response_kappa % get(p, this % xsData), TOL) @assertEqual(105.0_defReal, this % response_MT % get(p, this % xsData), TOL) p % isMG = .true. diff --git a/Tallies/TallyResponses/Tests/microResponse_test.f90 b/Tallies/TallyResponses/Tests/microResponse_test.f90 index ec1e810f9..462fac181 100644 --- a/Tallies/TallyResponses/Tests/microResponse_test.f90 +++ b/Tallies/TallyResponses/Tests/microResponse_test.f90 @@ -19,6 +19,7 @@ module microResponse_test type(microResponse) :: response_capture type(microResponse) :: response_fission type(microResponse) :: response_absorbtion + type(microResponse) :: response_kappa type(microResponse) :: response_MT type(testNeutronDatabase) :: xsData contains @@ -38,8 +39,8 @@ subroutine setUp(this) ! Allocate and initialise test nuclearData - ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission - call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal) + ! Cross-sections: Total eScattering IeScatter Capture Fission nuFission kappaFission + call this % xsData % build(6.0_defReal, 3.0_defReal, ZERO, 2.0_defReal, 1.0_defReal, 1.5_defReal, 9.0_defReal) ! Set dictionaries to initialise material call dictMat1 % init(1) @@ -97,6 +98,14 @@ subroutine setUp(this) call this % response_absorbtion % init(tempDict) call tempDict % kill() + ! Energy deposition + call tempDict % init(3) + call tempDict % store('type', 'microResponse') + call tempDict % store('MT', N_KAPPA) + call tempDict % store('material', 'Xenon') + call this % response_kappa % init(tempDict) + call tempDict % kill() + ! MT number call tempDict % init(2) call tempDict % store('type','microResponse') @@ -140,6 +149,7 @@ subroutine testGettingResponse(this) @assertEqual(0.5_defReal, this % response_fission % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_eScatter % get(p, this % xsData), TOL) @assertEqual(1.5_defReal, this % response_absorbtion % get(p, this % xsData), TOL) + @assertEqual(4.5_defReal, this % response_kappa % get(p, this % xsData), TOL) @assertEqual(8.0_defReal, this % response_MT % get(p, this % xsData), TOL) p % isMG = .true. diff --git a/Tallies/TallyResponses/macroResponse_class.f90 b/Tallies/TallyResponses/macroResponse_class.f90 index 44ba0c9f8..d5d488d84 100644 --- a/Tallies/TallyResponses/macroResponse_class.f90 +++ b/Tallies/TallyResponses/macroResponse_class.f90 @@ -115,6 +115,9 @@ subroutine build(self, MT) case(N_ABSORPTION) self % MT = macroAbsorbtion + case(N_KAPPA) + self % MT = macroKappaFission + case default self % mainData = .false. self % MT = MT diff --git a/Tallies/TallyResponses/microResponse_class.f90 b/Tallies/TallyResponses/microResponse_class.f90 index fd7c63808..9bd5bb72f 100644 --- a/Tallies/TallyResponses/microResponse_class.f90 +++ b/Tallies/TallyResponses/microResponse_class.f90 @@ -143,6 +143,9 @@ subroutine build(self, MT) case(N_ABSORPTION) self % MT = macroAbsorbtion + case(N_KAPPA) + self % MT = macroKappaFission + case default self % MT = MT self % mainData = .false. diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index f3d214b40..4cc8998c1 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -1501,8 +1501,9 @@ Example: :: * macroResponse: used to score macroscopic reaction rates - MT: MT number of the desired reaction. The options are: -1 (total), -2 (disappearance), - -3 (elastic scattering), -4 (total inelastic scattering), -6 (fission), -7 nu*fission), - -20 (total scattering), -21 (absorption), -22 (total non elastic, i.e., absorption + inelastic). + -3 (elastic scattering), -4 (total inelastic scattering), -6 (fission), -7 (nu*fission), + -20 (total scattering), -21 (absorption), -22 (total non elastic, i.e., absorption + inelastic), + -80 (kappa*fission). Additionally, all the MT numbers allowed by microResponse can be used here. Example: :: @@ -1516,7 +1517,7 @@ Example: :: * microResponse: used to score microscopic reaction rates - MT: MT number of the desired reaction. The options are: 1, 2, 3, 4, 5, 11, 16-25, 27-30, - 32-38, 41, 42, 44, 45, 51-90, 91, 101-109, 111-117, 203-207, 875-890. These MT numbers + 32-38, 41, 42, 44, 45, 51-90, 91, 101-109, 111-117, 203-207, 301, 875-890. These MT numbers are defined in the conventional way, i.e., following the ENDF standard - material: material name where to score the reaction. The material must be defined to include only one nuclide; its density could be anything, it doesn't