From 35c95030d56a1b997adda215ca1d13876390436a Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Wed, 17 Jan 2024 18:20:47 +0000 Subject: [PATCH 01/10] Adding majorant precalculation for DT --- .../aceDatabase/aceNeutronDatabase_class.f90 | 142 ++++++++++++++++-- .../baseMgNeutronDatabase_class.f90 | 48 ++++-- NuclearData/nuclearDatabase_inter.f90 | 18 ++- .../testNeutronDatabase_class.f90 | 15 ++ PhysicsPackages/eigenPhysicsPackage_class.f90 | 2 +- .../fixedSourcePhysicsPackage_class.f90 | 2 +- .../transportOperatorDT_class.f90 | 33 +++- .../transportOperatorFactory_func.f90 | 7 +- .../transportOperatorHT_class.f90 | 54 ++++--- .../transportOperatorST_class.f90 | 16 +- TransportOperator/transportOperator_inter.f90 | 8 +- 11 files changed, 282 insertions(+), 63 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 496dff341..af13bc421 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -3,7 +3,8 @@ module aceNeutronDatabase_class use numPrecision use endfConstants use universalVariables - use genericProcedures, only : fatalError, numToChar + use genericProcedures, only : fatalError, numToChar, concatenate, quickSort, & + removeDuplicatesSorted, binarySearch use dictionary_class, only : dictionary use RNG_class, only : RNG use charMap_class, only : charMap @@ -60,6 +61,8 @@ module aceNeutronDatabase_class !! 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 + !! eGrid -> 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 @@ -72,6 +75,8 @@ 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 :: eGrid real(defReal), dimension(2) :: Ebounds = ZERO integer(shortInt),dimension(:),allocatable :: activeMat @@ -81,16 +86,15 @@ module aceNeutronDatabase_class logical(defBool) :: hasDBRC = .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 + procedure :: initMajorant ! ceNeutronDatabase Procedures procedure :: energyBounds @@ -101,6 +105,10 @@ module aceNeutronDatabase_class procedure :: updateMicroXSs procedure :: getScattMicroMajXS + ! class Procedures + procedure :: init_urr + procedure :: init_DBRC + end type aceNeutronDatabase @@ -328,22 +336,26 @@ subroutine updateMajorantXS(self, E, rand) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand - integer(shortInt) :: i, matIdx + integer(shortInt) :: idx + 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) + idx = binarySearch(self % eGrid, E) - ! Update if needed - if( cache_materialCache(matIdx) % E_tot /= E) then - call self % updateTotalMatXS(E, matIdx, rand) - end if + if(idx <= 0) then + call fatalError(Here,'Failed to find energy: '//numToChar(E)//& + ' in unionised majorant grid') + end if + + associate(E_top => self % eGrid(idx + 1), E_low => self % eGrid(idx)) + f = (E - E_low) / (E_top - E_low) + end associate + + maj % xs = self % majorant(idx+1) * f + (ONE - f) * self % majorant(idx) - maj % xs = max(maj % xs, cache_materialCache(matIdx) % xss % total) - end do end associate end subroutine updateMajorantXS @@ -771,13 +783,111 @@ subroutine activate(self, 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 end subroutine activate + !! + !! Precomputes majorant cross section + !! + !! See nuclearDatabase documentation for details + !! + subroutine initMajorant(self, rand) + class(aceNeutronDatabase), intent(inout) :: self + class(RNG), intent(inout) :: rand + real(defReal), dimension(:), allocatable :: tmpGrid + integer(shortInt) :: i, j, matIdx, nNuc, nucIdx, isDone + type(intMap) :: nucMap + real(defReal) :: E, maj + integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 + + print '(A)', 'Building unionised energy grid' + + ! Initialise energy grid + matIdx = self % activeMat(1) + nucIdx = self % materials(matIdx) % nuclides(1) + tmpGrid = self % nuclides(nucIdx) % eGrid + + ! Add first nuclide to map to indicate it's been added + call nucMap % add(nucIdx, IN_SET) + + ! 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 grid + nucIdx = self % materials(matIdx) % nuclides(j) + isDone = nucMap % getOrDefault(nucIdx, NOT_PRESENT) + + ! If it's a new nuclide, add its energy grid to the unionised one + if (isDone /= IN_SET) then + + tmpGrid = concatenate(tmpGrid, self % nuclides(nucIdx) % eGrid) + + ! Sort and remove duplicates + call quickSort(tmpGrid) + tmpGrid = removeDuplicatesSorted(tmpGrid) + + ! Add nuclide to the set + call nucMap % add(nucIdx, IN_SET) + + end if + + end do + end do + + ! Save final grid + self % eGrid = tmpGrid + + print '(A)', 'Unionised energy grid has size: '//numToChar(size(self % eGrid))//& + '. Now building unionised majorant cross section' + + ! Allocate unionised majorant + allocate(self % majorant(size(self % eGrid))) + + ! Loop over all the energies + do i = 1, size(self % eGrid) + + ! Retrieve current energy + E = self % eGrid(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, rand) + maj = max(maj, cache_materialCache(matIdx) % xss % total) + + end do + + ! Save majorant for this energy + self % majorant(i) = maj + + end do + + print '(A)', 'Unionised majorant cross section completed' + + end subroutine initMajorant + !! !! Cast nuclearDatabase pointer to aceNeutronDatabase type pointer !! diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index 078c3f2b6..b7a19ac8f 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -7,6 +7,7 @@ module baseMgNeutronDatabase_class use charMap_class, only : charMap use dictionary_class, only : dictionary use dictParser_func, only : fileToDict + use RNG_class, only : RNG ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase @@ -70,6 +71,7 @@ module baseMgNeutronDatabase_class procedure :: kill procedure :: init procedure :: activate + procedure :: initMajorant ! Local interface procedure :: nGroups @@ -330,15 +332,41 @@ end subroutine init subroutine activate(self, activeMat) class(baseMgNeutronDatabase), intent(inout) :: self integer(shortInt), dimension(:), intent(in) :: activeMat + integer(shortInt) :: idx + + if(allocated(self % activeMats)) deallocate(self % activeMats) + self % activeMats = activeMat + + ! Initialies cross section cache + call cache_init(size(self % mats)) + + ! Store the material pointer in the material cache + !$omp parallel + do idx = 1,size(self % mats) + cache_materialCache(idx) % mat => self % mats(idx) + end do + !$omp end parallel + + end subroutine activate + + !! + !! Precomputes majorant cross section + !! + !! See nuclearDatabase documentation for details + !! + subroutine initMajorant(self, rand) + class(baseMgNeutronDatabase), intent(inout) :: self + class(RNG), intent(inout) :: rand integer(shortInt) :: g, i, idx real(defReal) :: xs integer(shortInt), parameter :: TOTAL_XS = 1 - if(allocated(self % activeMats)) deallocate(self % activeMats) - self % activeMats = activeMat + print '(A)', 'Building unionised majorant cross section' - ! Precalculate majorant xs for delta tracking + ! Allocate majorant allocate (self % majorant(self % nG)) + + ! Loop over energy groups do g = 1,self % nG xs = ZERO do i = 1,size(self % activeMats) @@ -348,17 +376,9 @@ subroutine activate(self, activeMat) self % majorant(g) = xs end do - ! Initialies cross section cache - call cache_init(size(self % mats)) + print '(A)', 'Unionised majorant cross section completed' - ! Store the material pointer in the material cache - !$omp parallel - do idx = 1,size(self % mats) - cache_materialCache(idx) % mat => self % mats(idx) - end do - !$omp end parallel - - end subroutine activate + end subroutine initMajorant !! !! Return number of energy groups in this database @@ -389,7 +409,7 @@ end function nGroups !! pure function baseMgNeutronDatabase_TptrCast(source) result(ptr) class(nuclearDatabase), pointer, intent(in) :: source - type(baseMgNeutronDatabase), pointer :: ptr + type(baseMgNeutronDatabase), pointer :: ptr select type(source) type is(baseMgNeutronDatabase) diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index 0c123df8d..c4a39352e 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -4,6 +4,7 @@ module nuclearDatabase_inter use dictionary_class, only : dictionary use particle_class, only : particle use charMap_class, only : charMap + use RNG_class, only : RNG ! Nuclear Data Handles use nuclideHandle_inter, only : nuclideHandle @@ -28,6 +29,7 @@ module nuclearDatabase_inter !! getMaterial -> returns a pointer to a material handle for given matIdx !! getNuclide -> returns a pointer to a nuclide handle for given nucIdx !! getReaction -> returns a pointer to a reaction for given matidx or nucIdx and MT number + !! initMajorant -> initialises the majorant cross section for delta tracking !! kill -> return to uninitialised state, clean memory !! type, public,abstract :: nuclearDatabase @@ -41,6 +43,7 @@ module nuclearDatabase_inter procedure(getMaterial), deferred :: getMaterial procedure(getNuclide), deferred :: getNuclide procedure(getReaction), deferred :: getReaction + procedure(initMajorant), deferred :: initMajorant procedure(kill), deferred :: kill end type nuclearDatabase @@ -196,7 +199,7 @@ end function getMaterial !! Allows to retrive an access to nuclide data for all databases types !! If database does not contain nuclides (e.g. MG data) just returns null() pointer !! - !! NOTE: This function can be used to inquire abou the presence of nucIdx in the database! + !! NOTE: This function can be used to inquire about the presence of nucIdx in the database! !! !! Args: !! nucIdx [in] -> nuclide index of required material @@ -225,7 +228,7 @@ end function getNuclide !! if MT < 0 then reaction is associated with material: idx -> matIdx !! if MT > 0 then reaction is associated with nuclide: idx -> nucIdx !! - !! NOTE: This function can be used to enquire abou the presence of data. If the data is + !! NOTE: This function can be used to enquire about the presence of data. If the data is !! not present null() pointer is always returned! !! !! Args: @@ -247,6 +250,17 @@ function getReaction(self, MT, idx) result(reac) class(reactionHandle),pointer :: reac end function getReaction + !! + !! Subroutine that precomputes the majorant cross section to be used by DT + !! + !! NOTE: Assumes that the nuclear database has been initialised and activated + !! + subroutine initMajorant(self, rand) + import :: nuclearDatabase, RNG + class(nuclearDatabase), intent(inout) :: self + class(RNG), intent(inout) :: rand + end subroutine initMajorant + !! !! Return to uninitialised state !! diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index 239f9384f..81f55b0b3 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -4,6 +4,7 @@ module testNeutronDatabase_class use particle_class, only : particle use dictionary_class, only : dictionary use charMap_class, only : charMap + use RNG_class, only : RNG ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase @@ -51,6 +52,7 @@ module testNeutronDatabase_class procedure :: getMaterial procedure :: getNuclide procedure :: getReaction + procedure :: initMajorant procedure :: kill ! Local Procedures @@ -261,6 +263,19 @@ function getReaction(self, MT, idx) result(reac) end function getReaction + !! + !! Initialise majorant for delta tracking + !! + !! See nuclearDatabase_inter for details + !! + subroutine initMajorant(self, rand) + class(testNeutronDatabase), intent(inout) :: self + class(RNG), intent(inout) :: rand + + ! Does nothing + + end subroutine initMajorant + !! !! Return to uninitialised state !! diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index aa89111c1..ff054afeb 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -490,7 +490,7 @@ subroutine init(self, dict) ! Build transport operator tempDict => dict % getDictPtr('transportOperator') - call new_transportOperator(self % transOp, tempDict) + call new_transportOperator(self % transOp, tempDict, self % particleType, self % pRNG) ! Initialise active & inactive tally Admins tempDict => dict % getDictPtr('inactiveTally') diff --git a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 index 35e2ff2d3..68ad1f8d8 100644 --- a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 +++ b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 @@ -406,7 +406,7 @@ subroutine init(self, dict) ! Build transport operator tempDict => dict % getDictPtr('transportOperator') - call new_transportOperator(self % transOp, tempDict) + call new_transportOperator(self % transOp, tempDict, self % particleType, self % pRNG) ! Initialise tally Admin tempDict => dict % getDictPtr('tally') diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 30d4b1e03..2c79b3959 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -9,10 +9,10 @@ module transportOperatorDT_class use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use rng_class, only : rng + use RNG_class, only : RNG ! Superclass - use transportOperator_inter, only : transportOperator + use transportOperator_inter, only : transportOperator, init_super => init ! Geometry interfaces use geometry_inter, only : geometry @@ -22,6 +22,7 @@ module transportOperatorDT_class use tallyAdmin_class, only : tallyAdmin ! Nuclear data interfaces + use nuclearDataReg_mod, only : ndReg_get => get use nuclearDatabase_inter, only : nuclearDatabase implicit none @@ -33,10 +34,15 @@ module transportOperatorDT_class type, public, extends(transportOperator) :: transportOperatorDT contains procedure :: transit => deltaTracking + ! Override procedure + procedure :: init end type transportOperatorDT contains + !! + !! Performs delta tracking until a real collision point is found + !! subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) class(transportOperatorDT), intent(inout) :: self class(particle), intent(inout) :: p @@ -91,7 +97,30 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) end do DTLoop call tally % reportTrans(p) + end subroutine deltaTracking + !! + !! Initialise DT transport operator + !! + !! See transportOperator_inter for more details + !! + subroutine init(self, dict, dataType, rand) + class(transportOperatorDT), intent(inout) :: self + class(dictionary), intent(in) :: dict + integer(shortInt), intent(in) :: dataType + class(RNG), intent(inout) :: rand + + ! Initialise superclass + call init_super(self, dict, dataType, rand) + + ! Get nuclear data pointer form the particle + self % xsData => ndReg_get(dataType) + + ! Precompute majorant cross section + call self % xsData % initMajorant(rand) + + end subroutine init + end module transportOperatorDT_class diff --git a/TransportOperator/transportOperatorFactory_func.f90 b/TransportOperator/transportOperatorFactory_func.f90 index a308ee01e..2fe492786 100644 --- a/TransportOperator/transportOperatorFactory_func.f90 +++ b/TransportOperator/transportOperatorFactory_func.f90 @@ -6,6 +6,7 @@ module transportOperatorFactory_func use numPrecision use genericProcedures, only : fatalError use dictionary_class, only : dictionary + use RNG_class, only : RNG ! Transport Operators use transportOperator_inter, only : transportOperator @@ -33,9 +34,11 @@ module transportOperatorFactory_func !! Allocate new allocatable transportOperator to a specific type !! If new is allocated it deallocates it !! - subroutine new_transportOperator(new, dict) + subroutine new_transportOperator(new, dict, dataType, rand) class(transportOperator),allocatable, intent(inout):: new class(dictionary), intent(in) :: dict + integer(shortInt), intent(in) :: dataType + class(RNG), intent(inout) :: rand character(nameLen) :: type character(100),parameter :: Here = 'new_transportOperator (transportOperatorFactory_func.f90)' @@ -62,7 +65,7 @@ subroutine new_transportOperator(new, dict) end select ! Initialise new transport operator - call new % init(dict) + call new % init(dict, dataType, rand) end subroutine new_transportOperator diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index f9cdfc378..0c4b31f8d 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -9,7 +9,7 @@ module transportOperatorHT_class use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use rng_class, only : rng + use RNG_class, only : RNG ! Tally interface use tallyCodes @@ -22,6 +22,7 @@ module transportOperatorHT_class use geometry_inter, only : geometry ! Nuclear data interfaces + use nuclearDataReg_mod, only : ndReg_get => get use nuclearDatabase_inter, only : nuclearDatabase implicit none @@ -31,29 +32,17 @@ module transportOperatorHT_class !! Transport operator that moves a particle with hybrid tracking !! type, public, extends(transportOperator) :: transportOperatorHT - !! Cutoff threshold between ST and DT - real(defReal) :: cutoff + real(defReal) :: cutoff ! Cutoff threshold between ST and DT contains - procedure :: init procedure :: transit => tracking_selection procedure, private :: deltaTracking procedure, private :: surfaceTracking + ! Override procedure + procedure :: init end type transportOperatorHT contains - subroutine init(self, dict) - class(transportOperatorHT), intent(inout) :: self - class(dictionary), intent(in) :: dict - - ! Initialise superclass - call init_super(self, dict) - - ! Initialise this class - call dict % getOrDefault(self % cutoff,'cutoff',0.9_defReal) - - end subroutine init - subroutine tracking_selection(self, p, tally, thisCycle, nextCycle) class(transportOperatorHT), intent(inout) :: self class(particle), intent(inout) :: p @@ -81,7 +70,9 @@ subroutine tracking_selection(self, p, tally, thisCycle, nextCycle) end subroutine tracking_selection - + !! + !! Performs delta tracking until a real collision point is found + !! subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) class(transportOperatorHT), intent(inout) :: self class(particle), intent(inout) :: p @@ -138,7 +129,9 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) call tally % reportTrans(p) end subroutine deltaTracking - + !! + !! Performs surface tracking until a collision point is found + !! subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) class(transportOperatorHT), intent(inout) :: self class(particle), intent(inout) :: p @@ -194,5 +187,30 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) end subroutine surfaceTracking + !! + !! Initialise HT operator from a dictionary + !! + !! See transportOperator_inter for more details + !! + subroutine init(self, dict, dataType, rand) + class(transportOperatorHT), intent(inout) :: self + class(dictionary), intent(in) :: dict + integer(shortInt), intent(in) :: dataType + class(RNG), intent(inout) :: rand + + ! Initialise superclass + call init_super(self, dict, dataType, rand) + + ! Retrieve DT-ST probability cutoff + call dict % getOrDefault(self % cutoff,'cutoff',0.9_defReal) + + ! Get nuclear data pointer form the particle + self % xsData => ndReg_get(dataType) + + ! Precompute majorant cross section for DT + call self % xsData % initMajorant(rand) + + end subroutine init + end module transportOperatorHT_class diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index c7eb24428..7c2927c58 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -1,5 +1,5 @@ !! -!! Transport operator for delta tracking +!! Transport operator for surface tracking !! module transportOperatorST_class use numPrecision @@ -28,7 +28,7 @@ module transportOperatorST_class private !! - !! Transport operator that moves a particle with delta tracking + !! Transport operator that moves a particle with surface tracking !! !! Sample Input Dictionary: !! trans { type transportOperatorST; cache 0;} @@ -37,13 +37,14 @@ module transportOperatorST_class logical(defBool) :: cache = .true. contains procedure :: transit => surfaceTracking + ! Override procedure procedure :: init end type transportOperatorST contains !! - !! Performs surface tracking + !! Performs surface tracking until a collision point is found !! subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) class(transportOperatorST), intent(inout) :: self @@ -108,13 +109,18 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) end subroutine surfaceTracking !! - !! Initialise surface operator from a dictionary + !! Initialise ST operator from a dictionary !! !! See transportOperator_inter for details !! - subroutine init(self, dict) + subroutine init(self, dict, dataType, rand) class(transportOperatorST), intent(inout) :: self class(dictionary), intent(in) :: dict + integer(shortInt), intent(in) :: dataType + class(RNG), intent(inout) :: rand + + ! Initialise superclass + call init_super(self, dict, dataType, rand) if (dict % isPresent('cache')) then call dict % get(self % cache, 'cache') diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index 8d0e30e68..53809254c 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -7,6 +7,8 @@ module transportOperator_inter use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary + use RNG_class, only : RNG + ! Geometry interfaces use geometryReg_mod, only : gr_geomPtr => geomPtr @@ -61,7 +63,7 @@ module transportOperator_inter end type transportOperator - ! Extandable procedures + ! Extendable procedures public :: init public :: kill @@ -120,9 +122,11 @@ end subroutine transport !! !! Initialise transport operator from dictionary and geometry !! - subroutine init(self, dict) + subroutine init(self, dict, dataType, rand) class(transportOperator), intent(inout) :: self class(dictionary), intent(in) :: dict + integer(shortInt), intent(in) :: dataType + class(RNG), intent(inout) :: rand ! Do nothing From 435aab702158dd3acc20c0d01d1c2c374e719be4 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Wed, 17 Jan 2024 18:45:07 +0000 Subject: [PATCH 02/10] Adjusting intTests adding majorant precalculation --- .../Tests/aceNeutronDatabase_iTest.f90 | 1 + .../aceDatabase/aceNeutronDatabase_class.f90 | 21 ++++++++++++++----- .../Tests/baseMgNeutronDatabase_iTest.f90 | 2 ++ .../baseMgNeutronDatabase_class.f90 | 15 ++++++++++--- NuclearData/nuclearDatabase_inter.f90 | 9 ++++---- .../testNeutronDatabase_class.f90 | 3 ++- 6 files changed, 38 insertions(+), 13 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 index a09f554d3..d1a8aef1a 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 @@ -76,6 +76,7 @@ subroutine test_aceNeutronDatabase() ptr => data call data % init(dataDict, ptr, silent = .true.) call data % activate([1,2]) + call data % initMajorant(p % pRNG, silent = .true.) !!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> !! Perform tests diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index af13bc421..966c5f6a8 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -795,16 +795,25 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, rand) + subroutine initMajorant(self, rand, silent) class(aceNeutronDatabase), intent(inout) :: self class(RNG), intent(inout) :: rand + logical(defBool), intent(in), optional :: silent + logical(defBool) :: loud real(defReal), dimension(:), allocatable :: tmpGrid integer(shortInt) :: i, j, matIdx, nNuc, nucIdx, isDone type(intMap) :: nucMap real(defReal) :: E, maj integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 - print '(A)', 'Building unionised energy grid' + ! Set build console output flag + if(present(silent)) then + loud = .not.silent + else + loud = .true. + end if + + if (loud) print '(A)', 'Building unionised energy grid' ! Initialise energy grid matIdx = self % activeMat(1) @@ -848,8 +857,10 @@ subroutine initMajorant(self, rand) ! Save final grid self % eGrid = tmpGrid - print '(A)', 'Unionised energy grid has size: '//numToChar(size(self % eGrid))//& - '. Now building unionised majorant cross section' + if (loud) then + print '(A)', 'Unionised energy grid has size: '//numToChar(size(self % eGrid))//& + '. Now building unionised majorant cross section' + end if ! Allocate unionised majorant allocate(self % majorant(size(self % eGrid))) @@ -884,7 +895,7 @@ subroutine initMajorant(self, rand) end do - print '(A)', 'Unionised majorant cross section completed' + if (loud) print '(A)', 'Unionised majorant cross section completed' end subroutine initMajorant diff --git a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 index 13aa797d0..1d8b017e8 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 @@ -75,6 +75,7 @@ subroutine testBaseMgNeutronDatabaseWithP0() call databaseDef % store('PN','P0') call database % init(databaseDef, data_ptr, silent = .true.) call database % activate([1]) + call database % initMajorant(p % pRNG, silent = .true.) ! Varify number of groups @assertEqual(4, database % nGroups()) @@ -201,6 +202,7 @@ subroutine testBaseMgNeutronDatabaseWithP1() call databaseDef % store('PN','P1') call database % init(databaseDef, data_ptr, silent = .true.) call database % activate([1]) + call database % initMajorant(p % pRNG, silent = .true.) ! Varify number of groups @assertEqual(4, database % nGroups()) diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index b7a19ac8f..e228b2388 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -354,14 +354,23 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, rand) + subroutine initMajorant(self, rand, silent) class(baseMgNeutronDatabase), intent(inout) :: self class(RNG), intent(inout) :: rand + logical(defBool), intent(in), optional :: silent + logical(defBool) :: loud integer(shortInt) :: g, i, idx real(defReal) :: xs integer(shortInt), parameter :: TOTAL_XS = 1 - print '(A)', 'Building unionised majorant cross section' + ! Set build console output flag + if(present(silent)) then + loud = .not.silent + else + loud = .true. + end if + + if (loud) print '(A)', 'Building unionised majorant cross section' ! Allocate majorant allocate (self % majorant(self % nG)) @@ -376,7 +385,7 @@ subroutine initMajorant(self, rand) self % majorant(g) = xs end do - print '(A)', 'Unionised majorant cross section completed' + if (loud) print '(A)', 'Unionised majorant cross section completed' end subroutine initMajorant diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index c4a39352e..6719ed8d6 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -255,10 +255,11 @@ end function getReaction !! !! NOTE: Assumes that the nuclear database has been initialised and activated !! - subroutine initMajorant(self, rand) - import :: nuclearDatabase, RNG - class(nuclearDatabase), intent(inout) :: self - class(RNG), intent(inout) :: rand + subroutine initMajorant(self, rand, silent) + import :: nuclearDatabase, RNG, defBool + class(nuclearDatabase), intent(inout) :: self + class(RNG), intent(inout) :: rand + logical(defBool), optional, intent(in) :: silent end subroutine initMajorant !! diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index 81f55b0b3..264c33e87 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -268,9 +268,10 @@ end function getReaction !! !! See nuclearDatabase_inter for details !! - subroutine initMajorant(self, rand) + subroutine initMajorant(self, rand, silent) class(testNeutronDatabase), intent(inout) :: self class(RNG), intent(inout) :: rand + logical(defBool), optional, intent(in) :: silent ! Does nothing From 18fbe55b11355df7ef2a8df4a248a3c9537021de Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Thu, 18 Jan 2024 16:52:14 +0000 Subject: [PATCH 03/10] Addressing comments --- .../Tests/aceNeutronDatabase_iTest.f90 | 3 +- .../aceDatabase/aceNeutronDatabase_class.f90 | 142 ++++++++++++------ .../ceNeutronData/ceNeutronDatabase_inter.f90 | 10 +- .../Tests/baseMgNeutronDatabase_iTest.f90 | 6 +- .../baseMgNeutronDatabase_class.f90 | 37 ++--- NuclearData/nuclearDatabase_inter.f90 | 20 +-- .../testNeutronDatabase_class.f90 | 18 +-- PhysicsPackages/eigenPhysicsPackage_class.f90 | 2 +- .../fixedSourcePhysicsPackage_class.f90 | 2 +- .../transportOperatorDT_class.f90 | 16 +- .../transportOperatorFactory_func.f90 | 11 +- .../transportOperatorHT_class.f90 | 18 +-- .../transportOperatorST_class.f90 | 9 +- TransportOperator/transportOperator_inter.f90 | 7 +- docs/User Manual.rst | 2 + 15 files changed, 151 insertions(+), 152 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 index d1a8aef1a..697cfcf9b 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 @@ -75,8 +75,7 @@ subroutine test_aceNeutronDatabase() ! Initialise data ptr => data call data % init(dataDict, ptr, silent = .true.) - call data % activate([1,2]) - call data % initMajorant(p % pRNG, silent = .true.) + call data % activate([1,2], silent = .true.) !!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> !! Perform tests diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 966c5f6a8..fd8bbceab 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -3,12 +3,13 @@ module aceNeutronDatabase_class use numPrecision use endfConstants use universalVariables - use genericProcedures, only : fatalError, numToChar, concatenate, quickSort, & - removeDuplicatesSorted, binarySearch - 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, concatenate, quickSort, & + 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 @@ -55,7 +56,8 @@ module aceNeutronDatabase_class !! Sample input: !! nuclearData { !! handles { - !! ce {type aceNeutronDatabase; DBRC (92238 94242); ures <1 or 0>; aceLibrary ;} } + !! ce { type aceNeutronDatabase; DBRC (92238 94242); ures < 1 or 0 >; + !! majorant < 1 or 0 >; aceLibrary ;} } !! !! Public Members: !! nuclides -> array of aceNeutronNuclides with data @@ -76,16 +78,18 @@ module aceNeutronDatabase_class type(aceNeutronNuclide),dimension(:),pointer :: nuclides => null() type(ceNeutronMaterial),dimension(:),pointer :: materials => null() real(defReal), dimension(:), allocatable :: majorant - real(defReal), dimension(:), allocatable :: eGrid + 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 + ! nuclearDatabase Procedures procedure :: kill procedure :: matNamesMap @@ -94,7 +98,6 @@ module aceNeutronDatabase_class procedure :: getReaction procedure :: init procedure :: activate - procedure :: initMajorant ! ceNeutronDatabase Procedures procedure :: energyBounds @@ -108,6 +111,7 @@ module aceNeutronDatabase_class ! class Procedures procedure :: init_urr procedure :: init_DBRC + procedure :: initMajorant end type aceNeutronDatabase @@ -297,7 +301,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 @@ -335,26 +339,46 @@ 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) :: idx + 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 - idx = binarySearch(self % eGrid, E) + ! Get majorant via the precomputed unionised cross section + if (self % hasMajorant) then + idx = binarySearch(self % eGridUnion, E) - if(idx <= 0) then - call fatalError(Here,'Failed to find energy: '//numToChar(E)//& - ' in unionised majorant grid') - end if + if(idx <= 0) then + call fatalError(Here,'Failed to find energy: '//numToChar(E)//& + ' in unionised majorant grid') + end if - associate(E_top => self % eGrid(idx + 1), E_low => self % eGrid(idx)) - f = (E - E_low) / (E_top - E_low) - end associate + 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 = self % majorant(idx+1) * f + (ONE - f) * self % majorant(idx) + 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 @@ -370,7 +394,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 @@ -410,7 +434,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), & @@ -445,7 +469,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) ) @@ -665,6 +689,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() @@ -773,9 +800,11 @@ 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) @@ -788,6 +817,36 @@ subroutine activate(self, activeMat) 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 !! @@ -795,25 +854,16 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, rand, silent) + subroutine initMajorant(self, loud) class(aceNeutronDatabase), intent(inout) :: self - class(RNG), intent(inout) :: rand - logical(defBool), intent(in), optional :: silent - logical(defBool) :: loud + logical(defBool), intent(in) :: loud real(defReal), dimension(:), allocatable :: tmpGrid integer(shortInt) :: i, j, matIdx, nNuc, nucIdx, isDone type(intMap) :: nucMap real(defReal) :: E, maj integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 - ! Set build console output flag - if(present(silent)) then - loud = .not.silent - else - loud = .true. - end if - - if (loud) print '(A)', 'Building unionised energy grid' + if (loud) print '(A)', 'Building CE unionised energy grid' ! Initialise energy grid matIdx = self % activeMat(1) @@ -855,21 +905,21 @@ subroutine initMajorant(self, rand, silent) end do ! Save final grid - self % eGrid = tmpGrid + self % eGridUnion = tmpGrid if (loud) then - print '(A)', 'Unionised energy grid has size: '//numToChar(size(self % eGrid))//& - '. Now building unionised majorant cross section' + print '(A)', 'CE unionised energy grid has size: '//numToChar(size(self % eGridUnion))//& + '. Now building CE unionised majorant cross section' end if ! Allocate unionised majorant - allocate(self % majorant(size(self % eGrid))) + allocate(self % majorant(size(self % eGridUnion))) ! Loop over all the energies - do i = 1, size(self % eGrid) + do i = 1, size(self % eGridUnion) ! Retrieve current energy - E = self % eGrid(i) + E = self % eGridUnion(i) ! Correct for energies higher or lower than the allowed boundaries if (E < self % EBounds(1)) E = self % EBounds(1) @@ -885,7 +935,7 @@ subroutine initMajorant(self, rand, silent) matIdx = self % activeMat(j) ! Calculate current material cross section and compare - call self % updateTotalMatXS(E, matIdx, rand) + call self % updateTotalMatXS(E, matIdx) maj = max(maj, cache_materialCache(matIdx) % xss % total) end do @@ -895,7 +945,7 @@ subroutine initMajorant(self, rand, silent) end do - if (loud) print '(A)', 'Unionised majorant cross section completed' + if (loud) print '(A)', 'CE unionised majorant cross section completed' end subroutine initMajorant diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index d0c8fd1d2..4e3168956 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -105,7 +105,7 @@ subroutine updateTotalMatXS(self, E, matIdx, rand) class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx - class(RNG), intent(inout) :: rand + class(RNG), optional, intent(inout) :: rand end subroutine updateTotalMatXS !! @@ -125,7 +125,7 @@ subroutine updateMajorantXS(self, E, rand) import :: ceNeutronDatabase, defReal, RNG class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E - class(RNG), intent(inout) :: rand + class(RNG), optional, intent(inout) :: rand end subroutine updateMajorantXS !! @@ -147,7 +147,7 @@ subroutine updateMacroXSs(self, E, matIdx, rand) class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx - class(RNG), intent(inout) :: rand + class(RNG), optional, intent(inout) :: rand end subroutine updateMacroXSs !! @@ -169,7 +169,7 @@ subroutine updateTotalXS(self, E, nucIdx, rand) class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: nucIdx - class(RNG), intent(inout) :: rand + class(RNG), optional, intent(inout) :: rand end subroutine updateTotalXS !! @@ -191,7 +191,7 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: nucIdx - class(RNG), intent(inout) :: rand + class(RNG), optional, intent(inout) :: rand end subroutine updateMicroXSs !! diff --git a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 index 1d8b017e8..21b1a5438 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/Tests/baseMgNeutronDatabase_iTest.f90 @@ -74,8 +74,7 @@ subroutine testBaseMgNeutronDatabaseWithP0() call databaseDef % init(1) call databaseDef % store('PN','P0') call database % init(databaseDef, data_ptr, silent = .true.) - call database % activate([1]) - call database % initMajorant(p % pRNG, silent = .true.) + call database % activate([1], silent = .true.) ! Varify number of groups @assertEqual(4, database % nGroups()) @@ -201,8 +200,7 @@ subroutine testBaseMgNeutronDatabaseWithP1() call databaseDef % init(1) call databaseDef % store('PN','P1') call database % init(databaseDef, data_ptr, silent = .true.) - call database % activate([1]) - call database % initMajorant(p % pRNG, silent = .true.) + call database % activate([1], silent = .true.) ! Varify number of groups @assertEqual(4, database % nGroups()) diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index e228b2388..0816c345d 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -2,12 +2,12 @@ module baseMgNeutronDatabase_class use numPrecision use endfConstants - use genericProcedures, only : fatalError, numToChar + use errors_mod, only : fatalError + use genericProcedures, only : numToChar use particle_class, only : particle use charMap_class, only : charMap use dictionary_class, only : dictionary use dictParser_func, only : fileToDict - use RNG_class, only : RNG ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase @@ -71,9 +71,9 @@ module baseMgNeutronDatabase_class procedure :: kill procedure :: init procedure :: activate - procedure :: initMajorant ! Local interface + procedure :: initMajorant procedure :: nGroups end type baseMgNeutronDatabase @@ -329,9 +329,11 @@ end subroutine init !! !! See nuclearDatabase documentation for details !! - subroutine activate(self, activeMat) + subroutine activate(self, activeMat, silent) class(baseMgNeutronDatabase), intent(inout) :: self integer(shortInt), dimension(:), intent(in) :: activeMat + logical(defBool), optional, intent(in) :: silent + logical(defBool) :: loud integer(shortInt) :: idx if(allocated(self % activeMats)) deallocate(self % activeMats) @@ -347,6 +349,16 @@ subroutine activate(self, activeMat) end do !$omp end parallel + ! Set build console output flag + if (present(silent)) then + loud = .not. silent + else + loud = .true. + end if + + ! Build unionised majorant + call self % initMajorant(loud) + end subroutine activate !! @@ -354,23 +366,14 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, rand, silent) + subroutine initMajorant(self, loud) class(baseMgNeutronDatabase), intent(inout) :: self - class(RNG), intent(inout) :: rand - logical(defBool), intent(in), optional :: silent - logical(defBool) :: loud + logical(defBool), intent(in) :: loud integer(shortInt) :: g, i, idx real(defReal) :: xs integer(shortInt), parameter :: TOTAL_XS = 1 - ! Set build console output flag - if(present(silent)) then - loud = .not.silent - else - loud = .true. - end if - - if (loud) print '(A)', 'Building unionised majorant cross section' + if (loud) print '(A)', 'Building MG unionised majorant cross section' ! Allocate majorant allocate (self % majorant(self % nG)) @@ -385,7 +388,7 @@ subroutine initMajorant(self, rand, silent) self % majorant(g) = xs end do - if (loud) print '(A)', 'Unionised majorant cross section completed' + if (loud) print '(A)', 'MG unionised majorant cross section completed' end subroutine initMajorant diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index 6719ed8d6..ff31af26f 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -4,7 +4,6 @@ module nuclearDatabase_inter use dictionary_class, only : dictionary use particle_class, only : particle use charMap_class, only : charMap - use RNG_class, only : RNG ! Nuclear Data Handles use nuclideHandle_inter, only : nuclideHandle @@ -29,7 +28,6 @@ module nuclearDatabase_inter !! getMaterial -> returns a pointer to a material handle for given matIdx !! getNuclide -> returns a pointer to a nuclide handle for given nucIdx !! getReaction -> returns a pointer to a reaction for given matidx or nucIdx and MT number - !! initMajorant -> initialises the majorant cross section for delta tracking !! kill -> return to uninitialised state, clean memory !! type, public,abstract :: nuclearDatabase @@ -43,7 +41,6 @@ module nuclearDatabase_inter procedure(getMaterial), deferred :: getMaterial procedure(getNuclide), deferred :: getNuclide procedure(getReaction), deferred :: getReaction - procedure(initMajorant), deferred :: initMajorant procedure(kill), deferred :: kill end type nuclearDatabase @@ -78,10 +75,11 @@ subroutine init(self, dict, ptr, silent) !! Errors: !! fatalError if activeMat contains materials not defined in the instance !! - subroutine activate(self, activeMat) - import :: nuclearDatabase, shortInt + subroutine activate(self, activeMat, silent) + import :: nuclearDatabase, shortInt, defBool class(nuclearDatabase), intent(inout) :: self integer(shortInt), dimension(:), intent(in) :: activeMat + logical(defBool), optional, intent(in) :: silent end subroutine activate !! @@ -250,18 +248,6 @@ function getReaction(self, MT, idx) result(reac) class(reactionHandle),pointer :: reac end function getReaction - !! - !! Subroutine that precomputes the majorant cross section to be used by DT - !! - !! NOTE: Assumes that the nuclear database has been initialised and activated - !! - subroutine initMajorant(self, rand, silent) - import :: nuclearDatabase, RNG, defBool - class(nuclearDatabase), intent(inout) :: self - class(RNG), intent(inout) :: rand - logical(defBool), optional, intent(in) :: silent - end subroutine initMajorant - !! !! Return to uninitialised state !! diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index 264c33e87..fa91ef9ed 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -52,7 +52,6 @@ module testNeutronDatabase_class procedure :: getMaterial procedure :: getNuclide procedure :: getReaction - procedure :: initMajorant procedure :: kill ! Local Procedures @@ -155,9 +154,10 @@ subroutine init(self, dict, ptr, silent) !! !! See nuclearDatabase_inter for details !! - subroutine activate(self, activeMat) + subroutine activate(self, activeMat, silent) class(testNeutronDatabase), intent(inout) :: self integer(shortInt), dimension(:), intent(in) :: activeMat + logical(defBool), optional, intent(in) :: silent ! Do nothing @@ -263,20 +263,6 @@ function getReaction(self, MT, idx) result(reac) end function getReaction - !! - !! Initialise majorant for delta tracking - !! - !! See nuclearDatabase_inter for details - !! - subroutine initMajorant(self, rand, silent) - class(testNeutronDatabase), intent(inout) :: self - class(RNG), intent(inout) :: rand - logical(defBool), optional, intent(in) :: silent - - ! Does nothing - - end subroutine initMajorant - !! !! Return to uninitialised state !! diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index ff054afeb..aa89111c1 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -490,7 +490,7 @@ subroutine init(self, dict) ! Build transport operator tempDict => dict % getDictPtr('transportOperator') - call new_transportOperator(self % transOp, tempDict, self % particleType, self % pRNG) + call new_transportOperator(self % transOp, tempDict) ! Initialise active & inactive tally Admins tempDict => dict % getDictPtr('inactiveTally') diff --git a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 index 68ad1f8d8..35e2ff2d3 100644 --- a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 +++ b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 @@ -406,7 +406,7 @@ subroutine init(self, dict) ! Build transport operator tempDict => dict % getDictPtr('transportOperator') - call new_transportOperator(self % transOp, tempDict, self % particleType, self % pRNG) + call new_transportOperator(self % transOp, tempDict) ! Initialise tally Admin tempDict => dict % getDictPtr('tally') diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 2c79b3959..3a1ab1609 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -5,11 +5,11 @@ module transportOperatorDT_class use numPrecision use universalVariables - use genericProcedures, only : fatalError, numToChar + use errors_mod, only : fatalError + use genericProcedures, only : numToChar use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use RNG_class, only : RNG ! Superclass use transportOperator_inter, only : transportOperator, init_super => init @@ -105,20 +105,12 @@ end subroutine deltaTracking !! !! See transportOperator_inter for more details !! - subroutine init(self, dict, dataType, rand) + subroutine init(self, dict) class(transportOperatorDT), intent(inout) :: self class(dictionary), intent(in) :: dict - integer(shortInt), intent(in) :: dataType - class(RNG), intent(inout) :: rand ! Initialise superclass - call init_super(self, dict, dataType, rand) - - ! Get nuclear data pointer form the particle - self % xsData => ndReg_get(dataType) - - ! Precompute majorant cross section - call self % xsData % initMajorant(rand) + call init_super(self, dict) end subroutine init diff --git a/TransportOperator/transportOperatorFactory_func.f90 b/TransportOperator/transportOperatorFactory_func.f90 index 2fe492786..549e2739b 100644 --- a/TransportOperator/transportOperatorFactory_func.f90 +++ b/TransportOperator/transportOperatorFactory_func.f90 @@ -4,9 +4,8 @@ module transportOperatorFactory_func use numPrecision - use genericProcedures, only : fatalError - use dictionary_class, only : dictionary - use RNG_class, only : RNG + use errors_mod, only : fatalError + use dictionary_class, only : dictionary ! Transport Operators use transportOperator_inter, only : transportOperator @@ -34,11 +33,9 @@ module transportOperatorFactory_func !! Allocate new allocatable transportOperator to a specific type !! If new is allocated it deallocates it !! - subroutine new_transportOperator(new, dict, dataType, rand) + subroutine new_transportOperator(new, dict) class(transportOperator),allocatable, intent(inout):: new class(dictionary), intent(in) :: dict - integer(shortInt), intent(in) :: dataType - class(RNG), intent(inout) :: rand character(nameLen) :: type character(100),parameter :: Here = 'new_transportOperator (transportOperatorFactory_func.f90)' @@ -65,7 +62,7 @@ subroutine new_transportOperator(new, dict, dataType, rand) end select ! Initialise new transport operator - call new % init(dict, dataType, rand) + call new % init(dict) end subroutine new_transportOperator diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index 0c4b31f8d..3ea03b9ae 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -5,11 +5,11 @@ module transportOperatorHT_class use numPrecision use universalVariables - use genericProcedures, only : fatalError, numToChar + use errors_mod, only : fatalError + use genericProcedures, only : numToChar use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use RNG_class, only : RNG ! Tally interface use tallyCodes @@ -89,7 +89,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) if (abs(majorant_inv) > huge(majorant_inv)) call fatalError(Here, "Majorant is 0") DTLoop:do - distance = -log( p% pRNG % get() ) * majorant_inv + distance = -log( p % pRNG % get() ) * majorant_inv ! Move partice in the geometry call self % geom % teleport(p % coords, distance) @@ -192,24 +192,16 @@ end subroutine surfaceTracking !! !! See transportOperator_inter for more details !! - subroutine init(self, dict, dataType, rand) + subroutine init(self, dict) class(transportOperatorHT), intent(inout) :: self class(dictionary), intent(in) :: dict - integer(shortInt), intent(in) :: dataType - class(RNG), intent(inout) :: rand ! Initialise superclass - call init_super(self, dict, dataType, rand) + call init_super(self, dict) ! Retrieve DT-ST probability cutoff call dict % getOrDefault(self % cutoff,'cutoff',0.9_defReal) - ! Get nuclear data pointer form the particle - self % xsData => ndReg_get(dataType) - - ! Precompute majorant cross section for DT - call self % xsData % initMajorant(rand) - end subroutine init diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index 7c2927c58..4b80605fc 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -5,11 +5,10 @@ module transportOperatorST_class use numPrecision use universalVariables - use genericProcedures, only : fatalError + use errors_mod, only : fatalError use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use RNG_class, only : RNG ! Superclass use transportOperator_inter, only : transportOperator, init_super => init @@ -113,14 +112,12 @@ end subroutine surfaceTracking !! !! See transportOperator_inter for details !! - subroutine init(self, dict, dataType, rand) + subroutine init(self, dict) class(transportOperatorST), intent(inout) :: self class(dictionary), intent(in) :: dict - integer(shortInt), intent(in) :: dataType - class(RNG), intent(inout) :: rand ! Initialise superclass - call init_super(self, dict, dataType, rand) + call init_super(self, dict) if (dict % isPresent('cache')) then call dict % get(self % cache, 'cache') diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index 53809254c..ddfdf5b6d 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -2,12 +2,11 @@ module transportOperator_inter use numPrecision use universalVariables - use genericProcedures, only : fatalError + use errors_mod, only : fatalError use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary - use RNG_class, only : RNG ! Geometry interfaces @@ -122,11 +121,9 @@ end subroutine transport !! !! Initialise transport operator from dictionary and geometry !! - subroutine init(self, dict, dataType, rand) + subroutine init(self, dict) class(transportOperator), intent(inout) :: self class(dictionary), intent(in) :: dict - integer(shortInt), intent(in) :: dataType - class(RNG), intent(inout) :: rand ! Do nothing diff --git a/docs/User Manual.rst b/docs/User Manual.rst index ef2841d4d..bbc12f9f4 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -725,6 +725,8 @@ from ACE files. resonance probability tables treatment * DBRC (*optional*, default = no DBRC): list of ZAIDs of nuclides for which DBRC has to be applied. +* majorant (*optional*, default = 0): 1 for true; 0 for false; flag to activate the + pre-construction of a unionised majorant cross section Example: :: From 66efce60878c244602210c068de81c892557b60c Mon Sep 17 00:00:00 2001 From: valeriaRaffuzzi <108435337+valeriaRaffuzzi@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:05:17 +0000 Subject: [PATCH 04/10] Update aceNeutronDatabase_class.f90 Add forgotten docs --- .../ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index fd8bbceab..8357b6a78 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -53,6 +53,9 @@ 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 { @@ -64,11 +67,12 @@ module aceNeutronDatabase_class !! materials -> array of ceNeutronMaterials with data !! Ebounds -> array with bottom (1) and top (2) energy bound !! majorant -> unionised majorant cross section - !! eGrid -> unionised energy grid + !! 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 From 8c82ffb4a05f165d9e041f1865951f9a7c9a75d0 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Wed, 24 Jan 2024 22:02:35 +0000 Subject: [PATCH 05/10] Speeding up unionised majorant cross section pre calculation --- .../aceDatabase/aceNeutronDatabase_class.f90 | 90 +++++++++++++------ .../baseMgNeutronDatabase_class.f90 | 4 +- 2 files changed, 65 insertions(+), 29 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index fd8bbceab..c475e4419 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -858,20 +858,15 @@ 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 - type(intMap) :: nucMap - real(defReal) :: E, maj + 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 - if (loud) print '(A)', 'Building CE unionised energy grid' - - ! Initialise energy grid - matIdx = self % activeMat(1) - nucIdx = self % materials(matIdx) % nuclides(1) - tmpGrid = self % nuclides(nucIdx) % eGrid - - ! Add first nuclide to map to indicate it's been added - call nucMap % add(nucIdx, IN_SET) + ! Find the size of the unionised energy grid (with duplicates) + ! Initialise size + sizeGrid = 0 ! Loop over active materials do i = 1, size(self % activeMat) @@ -883,33 +878,76 @@ subroutine initMajorant(self, loud) ! Loop over nuclides present in that material do j = 1, nNuc - ! Get index and check if it's already been added to the grid + ! Get index and check if it's already been added to the set nucIdx = self % materials(matIdx) % nuclides(j) - isDone = nucMap % getOrDefault(nucIdx, NOT_PRESENT) + isDone = nucSet % getOrDefault(nucIdx, NOT_PRESENT) - ! If it's a new nuclide, add its energy grid to the unionised one + ! If it's a new nuclide, add it to the set and find the size of its energy grid if (isDone /= IN_SET) then - tmpGrid = concatenate(tmpGrid, self % nuclides(nucIdx) % eGrid) + ! Add nuclide to the set + call nucSet % add(nucIdx, IN_SET) + + ! Update energy grid size + sizeGrid = sizeGrid + size(self % nuclides(nucIdx) % eGrid) - ! Sort and remove duplicates - call quickSort(tmpGrid) - tmpGrid = removeDuplicatesSorted(tmpGrid) + end if - ! Add nuclide to the set - call nucMap % add(nucIdx, IN_SET) + 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 - self % eGridUnion = tmpGrid + ! 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))//& - '. Now building CE unionised majorant cross section' + print '(A)', 'CE unionised energy grid has size: '//numToChar(size(self % eGridUnion)) end if ! Allocate unionised majorant @@ -945,7 +983,7 @@ subroutine initMajorant(self, loud) end do - if (loud) print '(A)', 'CE unionised majorant cross section completed' + if (loud) print '(A)', 'CE unionised majorant cross section calculation completed' end subroutine initMajorant diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index 0816c345d..a21280f5c 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -373,8 +373,6 @@ subroutine initMajorant(self, loud) real(defReal) :: xs integer(shortInt), parameter :: TOTAL_XS = 1 - if (loud) print '(A)', 'Building MG unionised majorant cross section' - ! Allocate majorant allocate (self % majorant(self % nG)) @@ -388,7 +386,7 @@ subroutine initMajorant(self, loud) self % majorant(g) = xs end do - if (loud) print '(A)', 'MG unionised majorant cross section completed' + if (loud) print '(A)', 'MG unionised majorant cross section calculation completed' end subroutine initMajorant From cb9a623144c24ce21f6de8448fe365e2df7bbcfc Mon Sep 17 00:00:00 2001 From: valeriaRaffuzzi <108435337+valeriaRaffuzzi@users.noreply.github.com> Date: Mon, 29 Jan 2024 16:52:52 +0000 Subject: [PATCH 06/10] Remove unused procedures in aceNeutronDatabase --- .../ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 37139a2f7..3fc8a3104 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -4,8 +4,7 @@ module aceNeutronDatabase_class use endfConstants use universalVariables use errors_mod, only : fatalError - use genericProcedures, only : numToChar, concatenate, quickSort, & - removeDuplicatesSorted, binarySearch + use genericProcedures, only : numToChar, removeDuplicatesSorted, binarySearch use dictionary_class, only : dictionary use RNG_class, only : RNG use charMap_class, only : charMap From 0072e914780f8cc57f0136a4003571759308faf7 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Wed, 31 Jan 2024 23:25:46 +0000 Subject: [PATCH 07/10] Typos and small fixes to DT majorant --- .../aceDatabase/aceNeutronDatabase_class.f90 | 22 +++++++++---------- NuclearData/nuclearDatabase_inter.f90 | 6 ++--- .../transportOperatorHT_class.f90 | 2 +- docs/User Manual.rst | 2 +- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 37139a2f7..21ac73d11 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -63,16 +63,16 @@ module aceNeutronDatabase_class !! majorant < 1 or 0 >; aceLibrary ;} } !! !! 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 + !! 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 + !! hasMajorant -> unionised majorant cross section flag !! !! Interface: !! nuclearData Interface @@ -694,7 +694,7 @@ subroutine init(self, dict, ptr, silent ) end do ! Read unionised majorant flag - call dict % getOrDefault(self % hasMajorant, 'majorant', .false.) + call dict % getOrDefault(self % hasMajorant, 'majorant', .true.) ! If on, initialise probability tables for ures if (self % hasUrr) then diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index ff31af26f..e3a8eeff2 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -171,7 +171,7 @@ end function matNamesMap !! !! Return pointer to material in a database !! - !! Allows to retrive an access to material data for all databases types + !! Allows to retrieve an access to material data for all databases types !! !! NOTE: This function can be used to inquire about the presence of matIdx in the database! !! @@ -194,7 +194,7 @@ end function getMaterial !! !! Return pointer to nuclide in a database !! - !! Allows to retrive an access to nuclide data for all databases types + !! Allows to retrieve an access to nuclide data for all databases types !! If database does not contain nuclides (e.g. MG data) just returns null() pointer !! !! NOTE: This function can be used to inquire about the presence of nucIdx in the database! @@ -218,7 +218,7 @@ end function getNuclide !! !! Return a pointer to a reaction !! - !! Allows to retrive an access to reaction data for all databases types + !! Allows to retrieve an access to reaction data for all databases types !! Reactions can be associated either with nuclides or materials. Thus, there is ambiguity !! whether material or nuclide should be asked to provide reaction data (using matIdx or nuIdx) !! diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index 3ea03b9ae..95ffaef94 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -91,7 +91,7 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) DTLoop:do distance = -log( p % pRNG % get() ) * majorant_inv - ! Move partice in the geometry + ! Move particle in the geometry call self % geom % teleport(p % coords, distance) ! If particle has leaked exit diff --git a/docs/User Manual.rst b/docs/User Manual.rst index bbc12f9f4..65da84db2 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -725,7 +725,7 @@ from ACE files. resonance probability tables treatment * DBRC (*optional*, default = no DBRC): list of ZAIDs of nuclides for which DBRC has to be applied. -* majorant (*optional*, default = 0): 1 for true; 0 for false; flag to activate the +* majorant (*optional*, default = 1): 1 for true; 0 for false; flag to activate the pre-construction of a unionised majorant cross section Example: :: From 04a37746f800d6cb2011526febefc0b56aa0d1b0 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Mon, 5 Feb 2024 17:17:57 +0000 Subject: [PATCH 08/10] Add routines to calculate majorant also with ures on --- .../Tests/aceNeutronDatabase_iTest.f90 | 9 + .../Tests/thermalScatteringData_iTest.f90 | 2 +- .../Tests/urrProbabilityTables_iTest.f90 | 2 +- .../aceDatabase/aceNeutronDatabase_class.f90 | 188 ++++++++++++++---- .../aceDatabase/aceNeutronNuclide_class.f90 | 35 ++-- .../urrProbabilityTables_class.f90 | 56 +++++- NuclearData/nuclearDatabase_inter.f90 | 1 + .../testNeutronDatabase_class.f90 | 1 - 8 files changed, 234 insertions(+), 60 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 index 697cfcf9b..69a08e08d 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 @@ -223,6 +223,15 @@ subroutine test_aceNeutronDatabase() p % E = 19.9_defReal @assertEqual(ONE, data % getMajorantXS(p)/0.21869599644_defReal , TOL) + ! Check that results are the same with on-the-fly majorant + data % hasMajorant = .false. + + p % E = 1.1E-6_defReal + @assertEqual(ONE, data % getMajorantXS(p) /4.4149556129495560_defReal , TOL) + + p % E = 19.9_defReal + @assertEqual(ONE, data % getMajorantXS(p)/0.21869599644_defReal , TOL) + !<><><><><><><><><><><><><><><><><><><><> ! Test getting Macroscopic XSs ! diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 index 447a5ed8d..61c75c77a 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 @@ -66,7 +66,7 @@ subroutine test_thermalScatteringData() ! 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 diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 index 149c6fa2d..6afb16e30 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/urrProbabilityTables_iTest.f90 @@ -57,7 +57,7 @@ subroutine test_urrProbabilityTables() ! Initialise data ptr => data call data % init(dataDict, ptr, silent = .true.) - call data % activate([1]) + call data % activate([1], silent = .true.) !!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> !! Perform tests diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 78280d7d0..a60c49b0a 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -52,9 +52,6 @@ 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 { @@ -112,8 +109,8 @@ module aceNeutronDatabase_class procedure :: getScattMicroMajXS ! class Procedures - procedure :: init_urr - procedure :: init_DBRC + procedure :: initUrr + procedure :: initDBRC procedure :: initMajorant end type aceNeutronDatabase @@ -446,16 +443,18 @@ subroutine updateTotalNucXS(self, E, nucIdx, rand) ! Check if the nuclide needs ures probability tables at this energy nucCache % needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2)) ! Check if the nuclide needs S(a,b) at this energy - nucCache % needsSabEl = (nuc % hasThData .and. E >= nuc % SabEl(1) .and. E < nuc % SabEl(2)) - nucCache % needsSabInel = (nuc % hasThData .and. E >= nuc % SabInel(1) .and. E < nuc % SabInel(2)) + nucCache % needsSabEl = (nuc % hasThData .and. E >= nuc % SabEl(1) .and. E <= nuc % SabEl(2)) + nucCache % needsSabInel = (nuc % hasThData .and. E >= nuc % SabInel(1) .and. E <= nuc % SabInel(2)) needsSab = (nucCache % needsSabEl .or. nucCache % needsSabInel) if (nucCache % needsUrr .or. needsSab) then call self % updateMicroXSs(E, nucIdx, rand) + else nucCache % E_tot = E call nuc % search(nucCache % idx, nucCache % f, E) nucCache % xss % total = nuc % totalXS(nucCache % idx, nucCache % f) + end if end associate @@ -489,17 +488,23 @@ subroutine updateMicroXSs(self, E, nucIdx, rand) ! Check if probability tables should be read if (nucCache % needsUrr) then associate(zaidCache => cache_zaidCache(self % nucToZaid(nucIdx))) + if (zaidCache % E /= E) then ! Save random number for temperature correlation zaidCache % xi = rand % get() zaidCache % E = E end if + call nuc % getUrrXSs(nucCache % xss, nucCache % idx, nucCache % f, E, zaidCache % xi) + end associate + elseif (nucCache % needsSabEl .or. nucCache % needsSabInel) then call nuc % getThXSs(nucCache % xss, nucCache % idx, nucCache % f, E) + else call nuc % microXSs(nucCache % xss, nucCache % idx, nucCache % f) + end if end associate @@ -642,11 +647,11 @@ subroutine init(self, dict, ptr, silent ) ! Initialise S(alpha,beta) tables if (idx /= 0 ) then call new_moderACE(ACE_Sab, name_file) - call self % nuclides(nucIdx) % init_Sab(ACE_Sab) + call self % nuclides(nucIdx) % initSab(ACE_Sab) end if ! Initialise probability tables - if (self % hasUrr) call self % nuclides(nucIdx) % init_urr(ACE) + if (self % hasUrr) call self % nuclides(nucIdx) % initUrr(ACE) ! Store nucIdx in the dictionary call nucSet % atSet(nucIdx, i) @@ -697,12 +702,12 @@ subroutine init(self, dict, ptr, silent ) ! If on, initialise probability tables for ures if (self % hasUrr) then - call self % init_urr() + call self % initUrr() end if ! If on, initialise DBRC if (self % hasDBRC) then - call self % init_DBRC(nucDBRC, nucSet, self % mapDBRCnuc) + call self % initDBRC(nucDBRC, nucSet, self % mapDBRCnuc) end if !! Clean up @@ -716,7 +721,7 @@ end subroutine init !! NOTE: compares the first 5 letters of the ZAID.TT. It would be wrong with isotopes !! with Z > 99 !! - subroutine init_urr(self) + subroutine initUrr(self) class(aceNeutronDatabase), intent(inout) :: self integer(shortInt) :: i, j character(nameLen) :: zaid @@ -747,12 +752,12 @@ subroutine init_urr(self) end if end do - end subroutine init_urr + end subroutine initUrr !! !! Checks through all nuclides, creates map with nuclides present and corresponding 0K nuclide !! - subroutine init_DBRC(self, nucDBRC, nucSet, map) + subroutine initDBRC(self, nucDBRC, nucSet, map) class(aceNeutronDatabase), intent(inout) :: self character(nameLen), dimension(:), intent(in) :: nucDBRC type(charMap), intent(in) :: nucSet @@ -796,7 +801,7 @@ subroutine init_DBRC(self, nucDBRC, nucSet, map) end do - end subroutine init_DBRC + end subroutine initDBRC !! !! Activate this nuclearDatabase @@ -830,23 +835,8 @@ subroutine activate(self, activeMat, silent) 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 + ! Precompute majorant cross section + call self % initMajorant(loud) end if @@ -861,11 +851,15 @@ 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 + integer(shortInt) :: i, j, k, matIdx, nNuc, nucIdx, isDone, & + sizeGrid, eIdx, nucIdxLast, eIdxLast, & + urrIdx type(intMap) :: nucSet - real(defReal) :: eRef, eNuc, E, maj + real(defReal) :: eRef, eNuc, E, maj, total, dens, urrMaj, & + nucXS, f, eMax, eMin + logical(defBool) :: needsUrr integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 + real(defReal), parameter :: NUDGE = 1.0e-06_defReal ! Find the size of the unionised energy grid (with duplicates) ! Initialise size @@ -894,6 +888,11 @@ subroutine initMajorant(self, loud) ! Update energy grid size sizeGrid = sizeGrid + size(self % nuclides(nucIdx) % eGrid) + ! If URR probability tables or S(a,b) tables are used, add their energy + ! boundary values to the grid to minimise interpolation errors + if (self % nuclides(nucIdx) % hasProbTab) sizeGrid = sizeGrid + 2 + if (self % nuclides(nucIdx) % hasThData) sizeGrid = sizeGrid + 3 + end if end do @@ -904,7 +903,8 @@ subroutine initMajorant(self, loud) tmpGrid = self % EBounds(2) ! Loop over the energy grid - do i = 1, sizeGrid + i = 1 + do while (i < sizeGrid) ! Loop over all nuclides in the set - here the value of the intMap is used as an energy index j = nucSet % begin() @@ -944,6 +944,76 @@ subroutine initMajorant(self, loud) ! Increment the energy index saved in the intMap for the nuclides whose energy was added call nucSet % add(nucIdxLast, eIdxLast + 1) + ! Loop over all nuclides again to add S(a,b) and ures energy boundaries to grid + j = nucSet % begin() + do while (j /= nucSet % end()) + + ! Retrieve energy in the grid and nuclide information + eMin = tmpGrid(i - 1) + eMax = tmpGrid(i) + nucIdx = nucSet % atKey(j) + + ! Check for URR probability tables + if (self % nuclides(nucIdx) % hasProbTab) then + + ! Lower energy boundary + E = self % nuclides(nucIdx) % urrE(1) + if (E >= eMin .and. E < eMax) then + tmpGrid(i) = E + tmpGrid(i + 1) = eMax + ! Update counter + i = i + 1 + end if + + ! Upper energy boundary + E = self % nuclides(nucIdx) % urrE(2) + if (E >= eMin .and. E < eMax) then + tmpGrid(i) = E + tmpGrid(i + 1) = eMax + ! Update counter + i = i + 1 + end if + + end if + + ! Check for Sab tables + if (self % nuclides(nucIdx) % hasThData) then + + ! Elastic upper energy boundary (NOTE: lower boundary is fixed) + E = self % nuclides(nucIdx) % SabEl(2) + if (E >= eMin .and. E < eMax ) then + tmpGrid(i) = E + tmpGrid(i + 1) = eMax + ! Update counter + i = i + 1 + end if + + ! Inelastic lower energy boundary + E = self % nuclides(nucIdx) % SabInel(1) + if (E >= eMin .and. E < eMax ) then + tmpGrid(i) = E + tmpGrid(i + 1) = eMax + ! Update counter + i = i + 1 + end if + + ! Inelastic upper energy boundary + E = self % nuclides(nucIdx) % SabInel(2) + if (E >= eMin .and. E < eMax ) then + tmpGrid(i) = E + tmpGrid(i + 1) = eMax + ! Update counter + i = i + 1 + end if + + end if + + j = nucSet % next(j) + + end do + + i = i + 1 + end do ! Save final grid and remove duplicates @@ -974,15 +1044,51 @@ subroutine initMajorant(self, loud) ! Get material index matIdx = self % activeMat(j) + total = ZERO + + ! Loop over nuclides + do k = 1, size(self % materials(matIdx) % nuclides) + dens = self % materials(matIdx) % dens(k) + nucIdx = self % materials(matIdx) % nuclides(k) + + associate (nuc => self % nuclides(nucIdx)) + + needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2)) + + ! Check if present nuclide uses URR tabes + if (needsUrr) then + + ! Find maximum URR table total XS + urrIdx = binarySearch(nuc % probTab % eGrid, E) + urrMaj = nuc % probTab % majorant(urrIdx) + + ! Check if URR tables contain xs or multiplicative factor + if (nuc % IFF == 1) then + call nuc % search(eIdx, f, E) + nucXS = nuc % totalXS(eIdx, f) * urrMaj + else + nucXS = urrMaj + end if + + else + call self % updateTotalNucXS(E, nucIdx) + nucXS = cache_nuclideCache(nucIdx) % xss % total + + end if + + end associate + + ! Update total material cross section + total = total + dens * nucXS + + end do - ! Calculate current material cross section and compare - call self % updateTotalMatXS(E, matIdx) - maj = max(maj, cache_materialCache(matIdx) % xss % total) + maj = max(maj, total) end do - ! Save majorant for this energy - self % majorant(i) = maj + ! Save majorant for this energy. Nudge it up to avoid small discrepancies + self % majorant(i) = maj * (ONE + NUDGE) end do diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 27481891b..2438812e3 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -103,9 +103,9 @@ module aceNeutronNuclide_class !! getThXSs -> return ceNeutronMicroXSs accounting for S(a,b) scattering treatment !! elScatteringMaj -> returns the elastic scattering majorant within an energy range given as input !! init -> build nuclide from aceCard - !! init_urr -> build list and mapping of nuclides to maintain temperature correlation + !! initUrr -> build list and mapping of nuclides to maintain temperature correlation !! when reading ures probability tables - !! init_Sab -> builds S(a,b) propertied from aceCard + !! initSab -> builds S(a,b) propertied from aceCard !! display -> print information about the nuclide to the console !! type, public, extends(ceNeutronNuclide) :: aceNeutronNuclide @@ -147,8 +147,8 @@ module aceNeutronNuclide_class procedure :: getThXSs procedure :: elScatteringMaj procedure :: init - procedure :: init_urr - procedure :: init_Sab + procedure :: initUrr + procedure :: initSab procedure :: display end type aceNeutronNuclide @@ -476,6 +476,7 @@ elemental subroutine getThXSs(self, xss, idx, f, E) ! Retrieve capture and fission cross sections as usual xss % capture = data(CAPTURE_XS, 2) * f + (ONE-f) * data(CAPTURE_XS, 1) + 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) @@ -741,14 +742,14 @@ subroutine init(self, ACE, nucIdx, database) ! Create a stack of MT reactions, devide them into ones that produce 2nd-ary ! particlues and pure absorbtion associate (MTs => ACE % getScatterMTs()) - do i=1,size(MTs) + do i = 1,size(MTs) if (MTs(i) == N_ANYTHING) cycle call scatterMT % push(MTs(i)) end do end associate associate (MTs => [ACE % getFissionMTs(), ACE % getCaptureMTs()]) - do i=1,size(MTs) + do i = 1,size(MTs) if(MTs(i) == N_FISSION) cycle ! MT=18 is already included with FIS block call absMT % push(MTs(i)) end do @@ -760,7 +761,7 @@ subroutine init(self, ACE, nucIdx, database) ! Load scattering reactions N = scatterMT % size() self % nMT = N - do i =1,N + do i = 1,N call scatterMT % pop(MT) self % MTdata(i) % MT = MT self % MTdata(i) % firstIdx = ACE % firstIdxMT(MT) @@ -783,7 +784,7 @@ subroutine init(self, ACE, nucIdx, database) end do ! Calculate Inelastic scattering XS - do i=1,self % nMT + do i = 1,self % nMT do j=1,size(self % mainData, 2) ! Find bottom and Top of the grid bottom = self % MTdata(i) % firstIdx @@ -804,7 +805,7 @@ subroutine init(self, ACE, nucIdx, database) self % mainData(TOTAL_XS, :) = sum(self % mainData(ESCATTER_XS:K,:),1) ! Load Map of MT -> local index of a reaction - do i=1,size(self % MTdata) + do i = 1,size(self % MTdata) call self % idxMT % add(self % MTdata(i) % MT, i) end do @@ -830,7 +831,7 @@ end subroutine init !! Args: !! ACE [inout] -> ACE card !! - subroutine init_urr(self, ACE) + subroutine initUrr(self, ACE) class(aceNeutronNuclide), intent(inout) :: self class(aceCard), intent(inout) :: ACE @@ -841,19 +842,24 @@ subroutine init_urr(self, ACE) ! Initialise probability tables call self % probTab % init(ACE) ! Check if probability tables were read correctly + if (allocated(self % probTab % eGrid)) then self % urrE = self % probTab % getEbounds() self % IFF = self % probTab % getIFF() + else ! Something went wrong! self % hasProbTab = .false. self % urrE = ZERO + end if + else self % urrE = ZERO + end if - end subroutine init_urr + end subroutine initUrr !! !! Initialise thermal scattering tables from ACE card @@ -865,14 +871,15 @@ end subroutine init_urr !! fatalError if the inelastic scattering S(a,b) energy grid starts at a !! lower energy than the nuclide energy grid !! - subroutine init_Sab(self, ACE) + subroutine initSab(self, ACE) class(aceNeutronNuclide), intent(inout) :: self class(aceSabCard), intent(inout) :: ACE - character(100), parameter :: Here = "init_Sab (aceNeutronNuclide_class.f90)" + character(100), parameter :: Here = "initSab (aceNeutronNuclide_class.f90)" ! Initialise S(a,b) class from ACE file call self % thData % init(ACE) self % hasThData = .true. + ! Initialise energy boundaries self % SabInel = self % thData % getEbounds('inelastic') self % SabEl = self % thData % getEbounds('elastic') @@ -882,7 +889,7 @@ subroutine init_Sab(self, ACE) call fatalError(Here, 'S(a,b) low energy boundary is lower than nuclide first energy point') end if - end subroutine init_Sab + end subroutine initSab !! !! A Procedure that displays information about the nuclide to the screen diff --git a/NuclearData/ceNeutronData/aceDatabase/urrProbabilityTables_class.f90 b/NuclearData/ceNeutronData/aceDatabase/urrProbabilityTables_class.f90 index 87f8ce93d..42e292aff 100644 --- a/NuclearData/ceNeutronData/aceDatabase/urrProbabilityTables_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/urrProbabilityTables_class.f90 @@ -35,8 +35,9 @@ module urrProbabilityTables_class !! ILF -> Inelatic competition flag (only used to check if inelatic scattering is zero) !! IOA -> Other absorption flag. NOTE: not used in the code !! IFF -> Multiplication factor flag - !! eGrid -> Energy grid - !! table -> Array of probability tables + !! majorant -> Array of maximum table total xs entry per each energy + !! eGrid -> Energy grid + !! table -> Array of probability tables !! !! type, public :: urrProbabilityTables @@ -46,12 +47,14 @@ module urrProbabilityTables_class integer(shortInt) :: ILF integer(shortInt) :: IOA integer(shortInt) :: IFF + real(defReal), dimension(:), allocatable :: majorant real(defReal), dimension(:), allocatable :: eGrid type(urrTable), dimension(:), allocatable :: table contains procedure :: init + procedure :: computeMajorant procedure :: kill procedure :: getEbounds procedure :: getIFF @@ -85,6 +88,9 @@ subroutine init(self, data) call fatalError(Here,'Probability tables cannot be build from '//data % myType()) end select + ! Find majorant + if (allocated(self % table)) call self % computeMajorant + end subroutine init !! @@ -215,21 +221,67 @@ subroutine buildFromACE(self, ACE) print '(A)', "Probability table discarded because CDF is not sorted" call self % kill() return + elseif (abs(self % table(i) % CDF(self % nTable) - ONE) > 1.0E-6_defReal) then print '(A)', "Probability table discarded because CDF does not end with 1.0 " call self % kill() return + elseif (self % table(i) % CDF(self % nTable) < ONE) then ! Adjust CDF if it is within tolerance but smaller than 1 due to floating point error self % table(i) % CDF(self % nTable) = ONE + elseif (any(self % table(i) % tot < ZERO) .or. any(self % table(i) % el < ZERO) .or. & any(self % table(i) % fiss < ZERO) .or. any(self % table(i) % capt < ZERO) ) then print '(A)', "Probability table discarded because negative cross-sections present " call self % kill() return + end if end do end subroutine buildFromACE + !! + !! Build majorant + !! + !! Finds the largest entry stored in each probability table (i.e., per each energy). + !! Each table entry is also compared with the equivalent entry for the following + !! energy; this is done so that interpolation between energies will not be needed + !! when using this majorant to compute the overall cross section majorant + !! + subroutine computeMajorant(self) + class(urrProbabilityTables), intent(inout) :: self + real(defReal) :: majorant + integer(shortInt) :: i, j + + ! Allocate majorant + allocate(self % majorant(self % nGrid)) + + ! Loop over energy grid + do i = 1,self % nGrid + + ! Initialise majorant value + majorant = ZERO + + ! Loop over table elements + do j = 1,self % nTable + + ! Update majorant if needed + majorant = max(majorant,self % table(i) % tot(j)) + + ! Check energy above (avoids the need for energy interpolation later) + if (i < self % nGrid) then + majorant = max(majorant,self % table(i + 1) % tot(j)) + end if + + end do + + ! Save majorant for this energy + self % majorant(i) = majorant + + end do + + end subroutine computeMajorant + end module urrProbabilityTables_class diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index e3a8eeff2..ba70e7bce 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -71,6 +71,7 @@ subroutine init(self, dict, ptr, silent) !! !! Args: !! activeMat [in] -> Array of matIdx of materials active in the simulation + !! silent [in] -> Optional. If set to .true. disables console output !! !! Errors: !! fatalError if activeMat contains materials not defined in the instance diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index fa91ef9ed..fcfacf7f2 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -4,7 +4,6 @@ module testNeutronDatabase_class use particle_class, only : particle use dictionary_class, only : dictionary use charMap_class, only : charMap - use RNG_class, only : RNG ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase From b97e613b935dc01983db5281725f8f135d929549 Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Mon, 5 Feb 2024 17:31:56 +0000 Subject: [PATCH 09/10] Typos and preventing out of bound array access --- .../ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 | 6 +++++- .../ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index a60c49b0a..a95dc0b44 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -949,7 +949,11 @@ subroutine initMajorant(self, loud) do while (j /= nucSet % end()) ! Retrieve energy in the grid and nuclide information - eMin = tmpGrid(i - 1) + if (i /= 1) then + eMin = tmpGrid(i - 1) + else + eMin = ZERO + end if eMax = tmpGrid(i) nucIdx = nucSet % atKey(j) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 2438812e3..95086e572 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -105,7 +105,7 @@ module aceNeutronNuclide_class !! init -> build nuclide from aceCard !! initUrr -> build list and mapping of nuclides to maintain temperature correlation !! when reading ures probability tables - !! initSab -> builds S(a,b) propertied from aceCard + !! initSab -> builds S(a,b) properties from aceCard !! display -> print information about the nuclide to the console !! type, public, extends(ceNeutronNuclide) :: aceNeutronNuclide From a65100108421aeaf131b499dba43af68da498c6e Mon Sep 17 00:00:00 2001 From: "V. Raffuzzi(vr339)" Date: Mon, 5 Feb 2024 19:01:33 +0000 Subject: [PATCH 10/10] Typos --- .../ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 | 2 +- .../ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 | 4 ++-- NuclearData/nuclearDatabase_inter.f90 | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index a95dc0b44..d1bc37e9d 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -1059,7 +1059,7 @@ subroutine initMajorant(self, loud) needsUrr = (nuc % hasProbTab .and. E >= nuc % urrE(1) .and. E <= nuc % urrE(2)) - ! Check if present nuclide uses URR tabes + ! Check if present nuclide uses URR tables if (needsUrr) then ! Find maximum URR table total XS diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 95086e572..343bc6b0e 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -739,8 +739,8 @@ subroutine init(self, ACE, nucIdx, database) ! Read data for MT reaction - ! Create a stack of MT reactions, devide them into ones that produce 2nd-ary - ! particlues and pure absorbtion + ! Create a stack of MT reactions, divide them into ones that produce 2nd-ary + ! particles and pure absorption associate (MTs => ACE % getScatterMTs()) do i = 1,size(MTs) if (MTs(i) == N_ANYTHING) cycle diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index ba70e7bce..0c8712fd0 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -221,7 +221,7 @@ end function getNuclide !! !! Allows to retrieve an access to reaction data for all databases types !! Reactions can be associated either with nuclides or materials. Thus, there is ambiguity - !! whether material or nuclide should be asked to provide reaction data (using matIdx or nuIdx) + !! whether material or nuclide should be asked to provide reaction data (using matIdx or nucIdx) !! !! This ambiguity is resolved by following convenction: !! if MT < 0 then reaction is associated with material: idx -> matIdx