DT majorant

NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90

@@ -75,7 +75,7 @@ subroutine test_aceNeutronDatabase()
     ! Initialise data
     ptr => data
     call data % init(dataDict, ptr, silent = .true.)
-    call data % activate([1,2])
+    call data % activate([1,2], silent = .true.)
 
     !!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
     !! Perform tests

NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90

@@ -3,11 +3,12 @@ module aceNeutronDatabase_class
   use numPrecision
   use endfConstants
   use universalVariables
-  use genericProcedures, only : fatalError, numToChar
-  use dictionary_class,  only : dictionary
-  use RNG_class,         only : RNG
-  use charMap_class,     only : charMap
-  use intMap_class,      only : intMap
+  use errors_mod,         only : fatalError
+  use genericProcedures,  only : numToChar, removeDuplicatesSorted, binarySearch
+  use dictionary_class,   only : dictionary
+  use RNG_class,          only : RNG
+  use charMap_class,      only : charMap
+  use intMap_class,       only : intMap
 
   ! Nuclear Data Interfaces
   use nuclearDatabase_inter,        only : nuclearDatabase
@@ -51,19 +52,26 @@ module aceNeutronDatabase_class
   !! It's possible to use probability tables in the unresolved resonance range if
   !! ures is included in the input file
   !!
+  !! NOTE: the unionised majorant is not calculated and used if probability tables 
+  !!       are on 
+  !!
   !! Sample input:
   !!   nuclearData {
   !!   handles {
-  !!   ce {type aceNeutronDatabase; DBRC (92238 94242); ures <1 or 0>; aceLibrary <nuclear data path> ;} }
+  !!   ce { type aceNeutronDatabase; DBRC (92238 94242); ures < 1 or 0 >;
+  !!        majorant < 1 or 0 >; aceLibrary <nuclear data path> ;} }
   !!
   !! Public Members:
   !!   nuclides   -> array of aceNeutronNuclides with data
   !!   materials  -> array of ceNeutronMaterials with data
   !!   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
   !!   nucToZaid  -> map to link nuclide index to zaid index
   !!   hasUrr     -> ures probability tables flag, it's false by default
   !!   hasDBRC    -> DBRC flag, it's false by default
+  !!   has majorant -> unionised majorant cross section flag
   !!
   !! Interface:
   !!   nuclearData Interface
@@ -72,24 +80,26 @@ module aceNeutronDatabase_class
   type, public, extends(ceNeutronDatabase) :: aceNeutronDatabase
     type(aceNeutronNuclide),dimension(:),pointer :: nuclides  => null()
     type(ceNeutronMaterial),dimension(:),pointer :: materials => null()
+    real(defReal), dimension(:), allocatable     :: majorant
+    real(defReal), dimension(:), allocatable     :: eGridUnion
     real(defReal), dimension(2)                  :: Ebounds   = ZERO
     integer(shortInt),dimension(:),allocatable   :: activeMat
 
     ! Probability tables data
     integer(shortInt),dimension(:),allocatable   :: nucToZaid
-    logical(defBool)                             :: hasUrr = .false.
+    logical(defBool)                             :: hasUrr  = .false.
     logical(defBool)                             :: hasDBRC = .false.
+    logical(defBool)                             :: hasMajorant = .false.
 
   contains
-    ! nuclearData Procedures
+
+    ! nuclearDatabase Procedures
     procedure :: kill
     procedure :: matNamesMap
     procedure :: getMaterial
     procedure :: getNuclide
     procedure :: getReaction
     procedure :: init
-    procedure :: init_urr
-    procedure :: init_DBRC
     procedure :: activate
 
     ! ceNeutronDatabase Procedures
@@ -101,6 +111,11 @@ module aceNeutronDatabase_class
     procedure :: updateMicroXSs
     procedure :: getScattMicroMajXS
 
+    ! class Procedures
+    procedure :: init_urr
+    procedure :: init_DBRC
+    procedure :: initMajorant
+
   end type aceNeutronDatabase
 
 
@@ -289,7 +304,7 @@ subroutine updateTotalMatXS(self, E, matIdx, rand)
     class(aceNeutronDatabase), intent(in) :: self
     real(defReal), intent(in)             :: E
     integer(shortInt), intent(in)         :: matIdx
-    class(RNG), intent(inout)             :: rand
+    class(RNG), optional, intent(inout)   :: rand
     integer(shortInt)                     :: i, nucIdx
     real(defReal)                         :: dens
 
@@ -327,23 +342,47 @@ end subroutine updateTotalMatXS
   subroutine updateMajorantXS(self, E, rand)
     class(aceNeutronDatabase), intent(in) :: self
     real(defReal), intent(in)             :: E
-    class(RNG), intent(inout)             :: rand
-    integer(shortInt)                     :: i, matIdx
+    class(RNG), optional, intent(inout)   :: rand
+    integer(shortInt)                     :: idx, i, matIdx
+    real(defReal)                         :: f
+    character(100), parameter :: Here = 'updateMajorantXS (aceNeutronDatabase_class.f90)'
 
     associate (maj => cache_majorantCache(1) )
       maj % E  = E
-      maj % xs = ZERO
 
-      do i = 1, size(self % activeMat)
-        matIdx = self % activeMat(i)
+      ! Get majorant via the precomputed unionised cross section
+      if (self % hasMajorant) then
+        idx = binarySearch(self % eGridUnion, E)
 
-        ! Update if needed
-        if( cache_materialCache(matIdx) % E_tot /= E) then
-          call self % updateTotalMatXS(E, matIdx, rand)
+        if(idx <= 0) then
+          call fatalError(Here,'Failed to find energy: '//numToChar(E)//&
+                               ' in unionised majorant grid')
         end if
 
-        maj % xs = max(maj % xs, cache_materialCache(matIdx) % xss % total)
-      end do
+        associate(E_top => self % eGridUnion(idx + 1), E_low  => self % eGridUnion(idx))
+          f = (E - E_low) / (E_top - E_low)
+        end associate
+
+        maj % xs = self % majorant(idx+1) * f + (ONE - f) * self % majorant(idx)
+
+      else ! Compute majorant on the fly
+
+        maj % xs = ZERO
+
+        ! Loop over materials
+        do i = 1, size(self % activeMat)
+          matIdx = self % activeMat(i)
+
+          ! Update if needed
+          if( cache_materialCache(matIdx) % E_tot /= E) then
+            call self % updateTotalMatXS(E, matIdx, rand)
+          end if
+
+          maj % xs = max(maj % xs, cache_materialCache(matIdx) % xss % total)
+        end do
+
+      end if
+
     end associate
 
   end subroutine updateMajorantXS
@@ -358,7 +397,7 @@ subroutine updateMacroXSs(self, E, matIdx, rand)
     class(aceNeutronDatabase), intent(in) :: self
     real(defReal), intent(in)             :: E
     integer(shortInt), intent(in)         :: matIdx
-    class(RNG), intent(inout)             :: rand
+    class(RNG), optional, intent(inout)   :: rand
     integer(shortInt)                     :: i, nucIdx
     real(defReal)                         :: dens
 
@@ -398,7 +437,7 @@ subroutine updateTotalNucXS(self, E, nucIdx, rand)
     class(aceNeutronDatabase), intent(in) :: self
     real(defReal), intent(in)             :: E
     integer(shortInt), intent(in)         :: nucIdx
-    class(RNG), intent(inout)             :: rand
+    class(RNG), optional, intent(inout)   :: rand
     logical(defBool)                      :: needsSab
 
     associate (nucCache => cache_nuclideCache(nucIdx), &
@@ -433,7 +472,7 @@ subroutine updateMicroXSs(self, E, nucIdx, rand)
     class(aceNeutronDatabase), intent(in) :: self
     real(defReal), intent(in)             :: E
     integer(shortInt), intent(in)         :: nucIdx
-    class(RNG), intent(inout)             :: rand
+    class(RNG), optional, intent(inout)   :: rand
 
     associate (nucCache => cache_nuclideCache(nucIdx), &
                nuc      => self % nuclides(nucIdx)     )
@@ -653,6 +692,9 @@ subroutine init(self, dict, ptr, silent )
       self % Ebounds(2) = min(self % Ebounds(2), self % nuclides(i) % eGrid(j))
     end do
 
+    ! Read unionised majorant flag
+    call dict % getOrDefault(self % hasMajorant, 'majorant', .false.)
+
     ! If on, initialise probability tables for ures
     if (self % hasUrr) then
        call self % init_urr()
@@ -761,23 +803,193 @@ end subroutine init_DBRC
   !!
   !! See nuclearDatabase documentation for details
   !!
-  subroutine activate(self, activeMat)
+  subroutine activate(self, activeMat, silent)
     class(aceNeutronDatabase), intent(inout)    :: self
     integer(shortInt), dimension(:), intent(in) :: activeMat
+    logical(defBool), optional, intent(in)      :: silent
+    logical(defBool)                            :: loud
 
     ! Load active materials
     if(allocated(self % activeMat)) deallocate(self % activeMat)
     self % activeMat = activeMat
 
     ! 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), 1, maxval(self % nucToZaid))
     else
-      call cache_init(size( self % materials), size(self % nuclides))
+      call cache_init(size(self % materials), size(self % nuclides))
+    end if
+
+    ! If unionised majorant cross section is requested, build it
+    if (self % hasMajorant) then
+
+      ! Set build console output flag
+      if (present(silent)) then
+        loud = .not. silent
+      else
+        loud = .true.
+      end if
+
+      ! Check if probability tables are on
+      if (self % hasUrr) then
+
+        ! Switch off majorant
+        self % hasMajorant = .false.
+
+        if (loud) then
+          print '(A)', 'Unionised majorant cross section will not be contructed &
+                        & due to the use of URR probability tables treatment'
+        end if
+
+      else
+
+        ! Precompute majorant cross section
+        call self % initMajorant(loud)
+
+      end if
+
     end if
 
   end subroutine activate
 
+  !!
+  !! Precomputes majorant cross section
+  !!
+  !! See nuclearDatabase documentation for details
+  !!
+  subroutine initMajorant(self, loud)
+    class(aceNeutronDatabase), intent(inout) :: self
+    logical(defBool), intent(in)             :: loud
+    real(defReal), dimension(:), allocatable :: tmpGrid
+    integer(shortInt)                        :: i, j, matIdx, nNuc, nucIdx, isDone, &
+                                                sizeGrid, eIdx, nucIdxLast, eIdxLast
+    type(intMap)                             :: nucSet
+    real(defReal)                            :: eRef, eNuc, E, maj
+    integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0
+
+    ! Find the size of the unionised energy grid (with duplicates)
+    ! Initialise size
+    sizeGrid = 0
+
+    ! Loop over active materials
+    do i = 1, size(self % activeMat)
+
+      ! Get current material index and number of nuclides in that material
+      matIdx = self % activeMat(i)
+      nNuc = size(self % materials(matIdx) % nuclides)
+
+      ! Loop over nuclides present in that material
+      do j = 1, nNuc
+
+        ! Get index and check if it's already been added to the set
+        nucIdx = self % materials(matIdx) % nuclides(j)
+        isDone = nucSet % getOrDefault(nucIdx, NOT_PRESENT)
+
+        ! If it's a new nuclide, add it to the set and find the size of its energy grid
+        if (isDone /= IN_SET) then
+
+          ! Add nuclide to the set
+          call nucSet % add(nucIdx, IN_SET)
+
+          ! Update energy grid size
+          sizeGrid = sizeGrid + size(self % nuclides(nucIdx) % eGrid)
+
+        end if
+
+      end do
+    end do
+
+    ! Allocate temporary grid vector and initialise to largest value allowed
+    allocate(tmpGrid(sizeGrid))
+    tmpGrid = self % EBounds(2)
+
+    ! Loop over the energy grid
+    do i = 1, sizeGrid
+
+      ! Loop over all nuclides in the set - here the value of the intMap is used as an energy index
+      j = nucSet % begin()
+      do while (j /= nucSet % end())
+
+        ! Retrieve energy in the grid and nuclide information
+        eRef    = tmpGrid(i)
+        nucIdx  = nucSet % atKey(j)
+        eIdx    = nucSet % atVal(j)
+
+        ! Check if we already added all the energy values for this nuclide
+        if (eIdx > size(self % nuclides(nucIdx) % eGrid)) then
+          j = nucSet % next(j)
+          cycle
+        end if
+
+        ! Get energy from nuclide grid
+        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
+          j = nucSet % next(j)
+          cycle
+        end if
+
+        ! Add energy value in the sorted grid, and save index of current nuclide
+        if (eNuc <= eRef) then
+          tmpGrid(i) = eNuc
+          nucIdxLast = nucIdx
+          eIdxLast   = eIdx
+        end if
+
+        j = nucSet % next(j)
+
+      end do
+
+      ! Increment the energy index saved in the intMap for the nuclides whose energy was added
+      call nucSet % add(nucIdxLast, eIdxLast + 1)
+
+    end do
+
+    ! Save final grid and remove duplicates
+    self % eGridUnion = removeDuplicatesSorted(tmpGrid)
+
+    if (loud) then
+      print '(A)', 'CE unionised energy grid has size: '//numToChar(size(self % eGridUnion))
+    end if
+
+    ! Allocate unionised majorant
+    allocate(self % majorant(size(self % eGridUnion)))
+
+    ! Loop over all the energies
+    do i = 1, size(self % eGridUnion)
+
+      ! Retrieve current energy
+      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)
+
+      ! Initialise majorant value for this energy
+      maj = ZERO
+
+      ! Loop over active materials
+      do j = 1, size(self % activeMat)
+
+        ! Get material index
+        matIdx = self % activeMat(j)
+
+        ! Calculate current material cross section and compare
+        call self % updateTotalMatXS(E, matIdx)
+        maj = max(maj, cache_materialCache(matIdx) % xss % total)
+
+      end do
+
+      ! Save majorant for this energy
+      self % majorant(i) = maj
+
+    end do
+
+    if (loud) print '(A)', 'CE unionised majorant cross section calculation completed'
+
+  end subroutine initMajorant
+
   !!
   !! Cast nuclearDatabase pointer to aceNeutronDatabase type pointer
   !!