Merge pull request CambridgeNuclear#1 from valeriaRaffuzzi/mathrad

Adding prompt and delayed nuFiss macroResponses

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

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

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 ,&

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