Fission energy release (#197)

Added fission energy release as a cross section

Author: ChasingNeutrons

IntegrationTestFiles/mgMat2

@@ -25,3 +25,5 @@ fission (1.0 0.0 0.0 0.0);
 chi (0.8 0.2 0.0 0.0);
 
 nu (2.3 0.0 0.0 0.0);
+
+kappa (202 193 180 200);

NuclearData/Reactions/reactionMG/Tests/fissionMG_test.f90

@@ -17,8 +17,7 @@ module fissionMG_test
   !!
   real(defReal),dimension(*),parameter :: nu = [2.3_defReal, 2.0_defReal, 1.3_defReal]
   real(defReal),dimension(*),parameter :: chi = [0.333333_defReal, 0.333333_defReal, 0.333334_defReal]
-
-
+  real(defReal),dimension(*),parameter :: kappa = [200.0_defReal, 203.0_defReal, 201.0_defReal]
 
 contains
 
@@ -63,6 +62,7 @@ subroutine fissionMG_Build_And_Functionality()
     ! Test Misc functionality
     @assertEqual(ZERO, reaction % releaseDelayed(1), TOL)
     @assertEqual(ZERO, reaction % sampleDelayRate(2, rand),TOL )
+    @assertEqual(202.27_defReal, reaction % getKappa(7),1E-5)
 
     ! Test Release
     @assertEqual(ZERO, reaction % releasePrompt(-2), TOL)
@@ -74,6 +74,19 @@ subroutine fissionMG_Build_And_Functionality()
     call dictT % kill()
     call reaction % kill()
 
+    ! Restart and test with kappa data included
+    call dictT % init(3)
+    call dictT % store('numberOfGroups',3)
+    call dictT % store('chi', chi)
+    call dictT % store('nu',nu)
+    call dictT % store('kappa', kappa)
+    
+    ! Build data Deck and initialise
+    data % dict => dictT
+    call reaction % init(data, macroFission)
+    @assertEqual(203.0_defReal, reaction % getKappa(2), 1E-5)
+
   end subroutine fissionMG_Build_And_Functionality
 
+
 end module fissionMG_test

NuclearData/Reactions/reactionMG/fissionMG_class.f90

@@ -22,18 +22,21 @@ module fissionMG_class
   !! A special type of MG reaction that contains data related to fission
   !!
   !! At the moment it has no information about delayed emissions, shall be introduced
-  !! at a latter date.
+  !! at a later date.
   !! 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, chi, or optionally kappa)]
+  !!   hasKappa -> flag for whether there is energy per fission data
   !!
   !! 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
+    logical(defBool)                         :: hasKappa = .false.
   contains
     ! Superclass procedures
     procedure :: init
@@ -46,13 +49,19 @@ module fissionMG_class
 
     ! Local procedures
     procedure :: buildFromDict
+    procedure :: getKappa
 
   end type fissionMG
 
   !!
   !! Reaction indices
   !!
-  integer(shortInt),parameter :: NU_DAT = 1, CHI_DAT = 2
+  integer(shortInt),parameter :: NU_DAT = 1, CHI_DAT = 2, KAPPA_DAT = 3
+  
+  !!
+  !! Default value of kappa used if no data is provided
+  !!
+  real(defReal), private, parameter :: KAPPA_DEFAULT = 202.27 ! [MeV]
 
 contains
 
@@ -216,9 +225,27 @@ subroutine buildFromDict(self, dict)
 
     ! Get number of groups
     call dict % get(nG, 'numberOfGroups')
+    
+    ! Check if energy per fission is present and allocate space appropriately
+    if (dict % isPresent('kappa')) then
+      
+      self % hasKappa = .true.
+      allocate(self % data(nG, 3))
+    
+      ! Get kappa
+      call dict % get(temp, 'kappa')
+      if(size(temp) /= ng) then
+        call fatalError(Here, 'Invalid number of values of kappa. Given: '// numToChar(size(temp)) // &
+                              ' Expected: ' // numToChar(nG))
+      end if
+      self % data(:,KAPPA_DAT) = temp
+
+    else
+      
+      self % hasKappa = .false.
+      allocate(self % data(nG, 2))
 
-    ! Allocate space
-    allocate(self % data(nG, 2))
+    end if
 
     ! Get nu
     call dict % get(temp, 'nu')
@@ -246,6 +273,22 @@ subroutine buildFromDict(self, dict)
 
   end subroutine buildFromDict
 
+  !!
+  !! Get the energy release from fission
+  !!
+  pure function getKappa(self, G) result(kappa)
+    class(fissionMG), intent(in)  :: self
+    integer(shortInt), intent(in) :: G
+    real(defReal)                 :: kappa
+
+    if (self % hasKappa) then
+      kappa = self % data(G, KAPPA_DAT)
+    else
+      kappa = KAPPA_DEFAULT
+    end if
+
+  end function getKappa
+
   !!
   !! Cast reactionHandle pointer to fissionMG pointer
   !!
@@ -270,5 +313,4 @@ pure function fissionMG_TptrCast(source) result(ptr)
 
   end function fissionMG_TptrCast
 
-
 end module fissionMG_class

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

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')
 

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

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
 

NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90

@@ -114,6 +114,7 @@ subroutine testBaseMgNeutronDatabaseWithP0()
     @assertEqual(1.0_defReal, xss % capture, TOL)
     @assertEqual(1.0_defReal, xss % fission, TOL)
     @assertEqual(2.3_defReal, xss % nuFission, TOL)
+    @assertEqual(202.0_defReal, xss % kappaXS, TOL)
 
     matClass => baseMgNeutronMaterial_CptrCast(database % getMaterial(1))
     @assertTrue(associated(matClass), "Type Ptr Cast has failed")
@@ -127,6 +128,7 @@ subroutine testBaseMgNeutronDatabaseWithP0()
     @assertEqual(4.0_defReal, xss % capture, TOL)
     @assertEqual(0.0_defReal, xss % fission, TOL)
     @assertEqual(0.0_defReal, xss % nuFission, TOL)
+    @assertEqual(0.0_defReal, xss % kappaXS, TOL)
 
     ! Get some invalid Materials
     mat => baseMgNeutronMaterial_TptrCast(database % getMaterial(0))
@@ -243,6 +245,7 @@ subroutine testBaseMgNeutronDatabaseWithP1()
     @assertEqual(1.0_defReal, xss % capture, TOL)
     @assertEqual(1.0_defReal, xss % fission, TOL)
     @assertEqual(2.3_defReal, xss % nuFission, TOL)
+    @assertEqual(202.0_defReal, xss % kappaXS, TOL)
 
     matClass => baseMgNeutronMaterial_CptrCast(database % getMaterial(1))
     @assertTrue(associated(matClass), "Type Ptr Cast has failed")
@@ -256,6 +259,7 @@ subroutine testBaseMgNeutronDatabaseWithP1()
     @assertEqual(4.0_defReal, xss % capture, TOL)
     @assertEqual(0.0_defReal, xss % fission, TOL)
     @assertEqual(0.0_defReal, xss % nuFission, TOL)
+    @assertEqual(0.0_defReal, xss % kappaXS, TOL)
 
     ! Get some invalid Materials
     mat => baseMgNeutronMaterial_TptrCast(database % getMaterial(0))