diff --git a/CMakeLists.txt b/CMakeLists.txt index 2861bb10d..b140c3b3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -135,6 +135,7 @@ add_subdirectory(UserInterface) add_subdirectory(PhysicsPackages) add_subdirectory(DataStructures) +add_subdirectory(Coupling) #################################################################################################### # Compile SCONE static library diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 index e90a69380..a3853a8f8 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 @@ -274,7 +274,7 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) do i = 1, n call fission % sampleOut(mu, phi, E_out, collDat % E, p % pRNG, lambda) - + ! Skip if a delayed particle is produced in prompt-only mode if (self % neglectDelayed .and. lambda < huge(lambda)) cycle @@ -420,7 +420,7 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) collDat % MT = self % nuc % invertInelastic(collDat % E, p % pRNG) reac => uncorrelatedReactionCE_CptrCast(self % xsData % getReaction(collDat % MT, collDat % nucIdx)) if (.not.associated(reac)) call fatalError(Here, "Failed to get scattering reaction") - + ! Scatter particle if (reac % inCMFrame()) then collDat % A = self % nuc % getMass() @@ -429,7 +429,7 @@ subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) call self % scatterInLAB(p, collDat, reac) end if - ! Apply weigth change using ingoing collision particle energy + ! Apply weight change using ingoing collision particle energy p % w = p % w * reac % release(collDat % E) end subroutine inelastic diff --git a/CollisionOperator/scatteringKernels_func.f90 b/CollisionOperator/scatteringKernels_func.f90 index c80df4afb..95952c93e 100644 --- a/CollisionOperator/scatteringKernels_func.f90 +++ b/CollisionOperator/scatteringKernels_func.f90 @@ -1,9 +1,10 @@ module scatteringKernels_func use numPrecision - use genericProcedures, only : rotateVector - use RNG_class, only : RNG - use particle_class, only : particle + use universalVariables, only : MINIMUM_ENERGY + use genericProcedures, only : rotateVector + use RNG_class, only : RNG + use particle_class, only : particle ! Nuclear Data use ceNeutronNuclide_inter, only : ceNeutronNuclide @@ -53,7 +54,8 @@ subroutine asymptoticScatter(E,mu,A) mu = (A*mu + 1)*sqrt(E_in/E)* inv_Ap1 ! Correct possible nuclear data shortcomings - if (mu > ONE) mu = ONE + mu = min(mu, ONE) + E = max(E, MINIMUM_ENERGY) end subroutine asymptoticScatter @@ -87,7 +89,8 @@ subroutine asymptoticInelasticScatter(E,mu,E_out,A) mu = mu * sqrt(E_out/E) + sqrt(E_in/E)* inv_Ap1 ! Correct possible nuclear data shortcomings - if (mu > ONE) mu = ONE + mu = min(mu, ONE) + E = max(E, MINIMUM_ENERGY) end subroutine asymptoticInelasticScatter diff --git a/Coupling/CMakeLists.txt b/Coupling/CMakeLists.txt new file mode 100644 index 000000000..044952249 --- /dev/null +++ b/Coupling/CMakeLists.txt @@ -0,0 +1,3 @@ +add_sources( ./couplingAdmin_class.f90) + +add_integration_tests( ./Tests/couplingAdmin_iTest.f90) diff --git a/Coupling/Tests/couplingAdmin_iTest.f90 b/Coupling/Tests/couplingAdmin_iTest.f90 new file mode 100644 index 000000000..d0a2f632b --- /dev/null +++ b/Coupling/Tests/couplingAdmin_iTest.f90 @@ -0,0 +1,180 @@ +module couplingAdmin_iTest + + use numPrecision + use universalVariables + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + use fieldFactory_func, only : new_field + use couplingAdmin_class, only : couplingAdmin + use funit + + implicit none + +contains + + !! + !! Coupling admin integration test + !! +@Test + subroutine test_coupling() + type(couplingAdmin) :: coupler + character(pathLen), parameter :: path = './IntegrationTestFiles/Coupling/test_coupling' + character(pathLen), parameter :: sendPath = './IntegrationTestFiles/Coupling/toDriver.dat' + character(pathLen), parameter :: recvPath = './IntegrationTestFiles/Coupling/toSCONE.dat' + character(pathLen), parameter :: tallyPath = './IntegrationTestFiles/Coupling/sconeOut.m' + type(dictionary) :: dict + integer(shortInt) :: it, unit + logical(defBool) :: exists, correctFile + character(nameLen) :: line + + ! Create a recv file to prevent SCONE from hanging + open(newunit=unit, file=recvPath, status="replace", action="write") + write(unit, '(A)') 'SIGUSR1' + close(unit) + + ! Attempt a coupling update before initialisation: no files should be produced + ! or deleted + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Initialise + call fileToDict(dict, path) + call coupler % init(dict) + + ! Feed in different iteration numbers and see whether files are produced + ! This one should not result in any information exchange + it = 1 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! This one should be coupled + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGUSR1 + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGUSR1') + @assertTrue(correctFile) + + ! Reset the files + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + + ! This one should end coupling due to hitting the maximum iteration + ! No recv file is necessary + it = 10 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGTERM + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGTERM') + @assertTrue(correctFile) + + ! Reset the files + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=recvPath, status="replace", action="write") + write(unit, '(A)') 'SIGTERM' + close(unit) + + ! Coupling should not occur + it = 15 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Reinitialise the coupler, ensuring that no coupling occurs successfully + ! until reinitialisation + call coupler % kill() + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertTrue(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + call coupler % init(dict) + + ! Given the recv file, coupling should occur on this iteration, but not on the 10th + it = 5 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertTrue(exists) + inquire(file=sendPath, exist=exists) + @assertTrue(exists) + + ! Send file should contain SIGUSR1 + open(newunit=unit, file=sendPath, status='old', action='read') + read(unit,'(A)') line + close(unit) + correctFile = (trim(line) == 'SIGUSR1') + @assertTrue(correctFile) + + ! Reset the files - no need for the recv file + open(newunit=unit, file=tallyPath, status='old') + close(unit, status='delete') + open(newunit=unit, file=sendPath, status='old') + close(unit, status='delete') + + it = 10 + call coupler % couple(it) + + inquire(file=recvPath, exist=exists) + @assertFalse(exists) + inquire(file=tallyPath, exist=exists) + @assertFalse(exists) + inquire(file=sendPath, exist=exists) + @assertFalse(exists) + + ! Kill coupling admin + call coupler % kill() + + end subroutine test_coupling + +end module couplingAdmin_iTest diff --git a/Coupling/couplingAdmin_class.f90 b/Coupling/couplingAdmin_class.f90 new file mode 100644 index 000000000..012319051 --- /dev/null +++ b/Coupling/couplingAdmin_class.f90 @@ -0,0 +1,468 @@ +module couplingAdmin_class + + use numPrecision + use universalVariables, only : nameDensity, nameTemperature, NO_TEMPERATURE + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + use outputFile_class, only : outputFile + use errors_mod, only : fatalError + use display_func, only : statusMsg + use genericProcedures, only : numToChar + use timer_mod, only : c_usleep + + use field_inter, only : field + use fieldFactory_func, only : new_field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast + use tallyAdmin_class, only : tallyAdmin + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_hasField => hasField + + implicit none + private + + character(nameLen), parameter, dimension(2) :: ALLOWABLE_FIELDS = ['temperature',& + 'density '] + character(nameLen), parameter :: END_SIGNAL = "SIGTERM", & + CONTINUE_SIGNAL = "SIGUSR1" + + !! + !! Object to manage coupling between radiation and external physics solvers/driver scripts + !! + !! This object is responsible for managing the status of coupled calculations. It should do the following: + !! - say whether coupling is happening (default: no) + !! - know the means of communication by which coupling occurs (currently files are used as flags) + !! - send a signal to the driver that it has finished a number of iterations + !! - wait for the driver/other physics to finish their calculations + !! - restart on receiving a signal from the driver + !! - create any tallies necessary for coupling (such as power tallies) + !! - output the results of those tallies to an appropriate location + !! - know the location of fields provided by other physics and update them in SCONE + !! - possibly provide checks on validity of information fed back to SCONE + !! - decide whether to continue coupling during active cycles + !! + !! Private members: + !! isCoupled -> flag for whether coupling is being performed + !! endWhenActive -> flag to end coupling when moving to active cycles (false by default) + !! commFileIn -> path to the file which triggers a SCONE calculation by the driver + !! commFileOut -> path to the file which SCONE outputs to the driver on completing a calculation + !! outputFile -> path to the file which SCONE outputs containing results to be used by the driver + !! updateFreq -> frequency with which SCONE pauses to exchange data with the driver + !! maxIt -> maximum number of iterations before finishing coupling + !! maxTemperature -> Assumed maximum temperature present in the problem + !! tally -> Tally admin containing necessary tallies for the driver + !! fieldPaths -> paths to files which SCONE reads to update the given fields + !! fieldNames -> names of fields which are updated by other physics + !! + !! Public interface: + !! init -> initialises the coupling setup from a dictionary + !! kill -> clean up + !! couple -> performs the main coupling logic + !! endCoupling -> deactivates coupling + !! moveToActive -> informs admin of shift to active cycles, ending coupling if needs be + !! + !! doCoupling -> returns a logical if coupling is happening + !! attachTally -> attach coupling tallies to another tallyAdmin used by a physics package + !! + !! Private procedures: + !! doUpdate -> returns a logical if a physics update should be performed on a given iteration. + !! maxTemperatureCheck -> checks whether the maximum field temperature is lower than that assumed initially + !! + !! outputTallies -> outputs tally data to be read by other physics + !! updateFields -> updates the relevant physics fields before resuming a calculation + !! + !! waitForSignal -> pauses the calculation until a 'resume' signal is received + !! signalIterationOver -> signals that an iteration is complete + !! signalCalculationOver -> signals that the coupled calculation is complete + !! + type, public :: couplingAdmin + private + logical(defBool) :: isCoupled = .false. + logical(defBool) :: endWhenActive = .false. + character(pathLen) :: commFileIn = '' + character(pathLen) :: commFileOut = '' + character(pathLen) :: outputFile = '' + character(nameLen) :: outputFormat = '' + integer(shortInt) :: updateFreq = 0 + integer(shortInt) :: maxIt = huge(shortInt) + real(defReal) :: maxTemperature = NO_TEMPERATURE + type(tallyAdmin), pointer :: tally + character(pathLen), dimension(:), allocatable :: fieldPaths + character(nameLen), dimension(:), allocatable :: fieldNames + contains + procedure :: init + procedure :: kill + procedure :: couple + procedure :: endCoupling + procedure :: moveToActive + + ! Status checks + procedure :: doCoupling + procedure, private :: checkMaxTemperature + procedure, private :: doUpdate + + ! Results handling and manipulation + procedure :: attachTally + procedure, private :: outputTallies + procedure, private :: updateFields + + ! Communication + procedure, private :: waitForSignal + procedure, private :: signalIterationOver + procedure, private :: signalCalculationOver + + end type couplingAdmin + +contains + + !! + !! Initialises coupling from a dictionary + !! + subroutine init(self, dict) + class(couplingAdmin), intent(inout) :: self + type(dictionary), intent(in) :: dict + integer(shortInt) :: i, nFields, j + logical(defBool) :: found + class(dictionary), pointer :: tallyDict + type(outputFile) :: test_out + character(100), parameter :: Here ='init (couplingAdmin_class.f90)' + + self % isCoupled = .true. + + ! Read communication details + call dict % get(self % commFileIn, 'receiveFile') + call dict % get(self % commFileOut, 'sendFile') + + ! Read whether to end when active + call dict % getOrDefault(self % endWhenActive, 'endActive', .false.) + + ! Read how often to exchange data + call dict % get(self % updateFreq, 'updateFreq') + if (self % updateFreq < 1) call fatalError(Here,& + 'Physics update frequency is silly: '//numToChar(self % updateFreq)) + + ! Read maximum number of iterations over which to perform coupling + call dict % getOrDefault(self % maxIt, 'maxIt', huge(shortInt)) + if (self % maxIt < 1) call fatalError(Here, & + 'Maximum iteration number is silly: '//numToChar(self % maxIt)) + + ! Read tally info needed by other solvers + ! Should this be modified to ensure MPI sync? + tallyDict => dict % getDictPtr('tally') + allocate(self % tally) + call self % tally % init(tallyDict) + + ! Read maximum temperature that may occur + call dict % getOrDefault(self % maxTemperature, 'maxTemp', NO_TEMPERATURE) + + ! Read the output file name + call dict % get(self % outputFile, 'outputFile') + + ! Initialise output file before calculation (so mistake in format will be caught early) + call dict % getOrDefault(self % outputFormat, 'outputFormat', 'asciiMATLAB') + call test_out % init(self % outputFormat) + + ! Read the field update location(s) + call dict % get(self % fieldPaths, 'fieldPaths') + + ! Read the field name(s) + call dict % get(self % fieldNames, 'fieldNames') + + ! Ensure the paths correspond to the name + if (size(self % fieldPaths) /= size(self % fieldNames)) call fatalError(Here,& + 'The number of field paths must match the number of field names.') + + ! Ensure the field names are unique and allowable + nFields = size(self % fieldNames) + do i = 1, nFields + + found = .false. + do j = 1, size(ALLOWABLE_FIELDS) + found = found .or. (trim(self % fieldNames(i)) == trim(ALLOWABLE_FIELDS(j))) + end do + + if (.not. found) then + print '(A)', "ALLOWABLE FIELDS:" + print '(A)', ALLOWABLE_FIELDS + call fatalError(Here,'Field name is invalid. See list above.') + end if + + do j = i+1, nFields + if (trim(self % fieldNames(i)) == trim(self % fieldNames(j))) then + call fatalError(Here,'Field names must be unique: '//trim(self % fieldNames(i))) + end if + end do + + end do + + end subroutine init + + !! + !! Clean-up + !! + subroutine kill(self) + class(couplingAdmin), intent(inout) :: self + + self % isCoupled = .false. + self % endWhenActive = .false. + self % updateFreq = 0 + self % commFileIn = '' + self % commFileOut = '' + self % outputFile = '' + self % outputFormat = '' + self % maxIt = huge(shortInt) + self % maxTemperature = NO_TEMPERATURE + call self % tally % kill() + if (allocated(self % fieldPaths)) deallocate(self % fieldPaths) + if (allocated(self % fieldNames)) deallocate(self % fieldNames) + + end subroutine kill + + !! + !! Subroutine which performs coupling logic during transport. + !! Evaluates whether a physics update should be performed. + !! If so, outputs tally information, cleans up tallies, signals and waits, + !! and updates fields. Deactivates coupling if later than the maximum + !! iteration. + !! + subroutine couple(self, it) + class(couplingAdmin), intent(inout) :: self + integer(shortInt), intent(in) :: it + + if (.not. self % doCoupling()) return + if (.not. self % doUpdate(it)) return + + call self % outputTallies() + call self % tally % resetMemory() + + if (it < self % maxIt) then + call self % signalIterationOver() + call self % waitForSignal() + call self % updateFields() + else + call self % signalCalculationOver() + self % isCoupled = .false. + end if + + end subroutine couple + + !! + !! Compares the maximum temperature field value with the assumed maximum + !! temperature given in the input. If this is exceeded, emits a warning + !! concerning the validity of the majorant and the resulting simulation + !! outputs. + !! + subroutine checkMaxTemperature(self, temp) + class(couplingAdmin), intent(in) :: self + real(defReal), intent(in) :: temp + + if (temp > self % maxTemperature) then + call statusMsg('WARNING: Coupled field temperature exceeds the assumed maximum temperature '& + //'for the simulation. If using delta tracking, results may be incorrect! '& + //'Assumed maximum: '//numToChar(self % maxTemperature)//'K. '& + //'Actual maximum: '//numToChar(temp)//'K.') + end if + + end subroutine checkMaxTemperature + + !! + !! Returns whether coupling should be occurring or not + !! + function doCoupling(self) result(isCoupled) + class(couplingAdmin), intent(in) :: self + logical(defBool) :: isCoupled + + isCoupled = self % isCoupled + + end function doCoupling + + !! + !! End coupling + !! + subroutine endCoupling(self) + class(couplingAdmin), intent(inout) :: self + + if (self % isCoupled) then + self % isCoupled = .false. + call self % signalCalculationOver() + end if + + end subroutine endCoupling + + !! + !! Signals to admin that active cycles have begun + !! Decides whether coupling should end. + !! If it continues, attach the active tally + !! + subroutine moveToActive(self, tallyPtr) + class(couplingAdmin), intent(inout) :: self + type(tallyAdmin), pointer, intent(inout) :: tallyPtr + + if (self % endWhenActive) then + call self % endCoupling() + else + call self % attachTally(tallyPtr) + end if + + end subroutine moveToActive + + !! + !! Returns whether to perform a physics update based on iteration number + !! and update frequency + !! + function doUpdate(self, it) result(update) + class(couplingAdmin), intent(in) :: self + integer(shortInt), intent(in) :: it + logical(defBool) :: update + + update = mod(it, self % updateFreq) == 0 + + end function doUpdate + + !! + !! Update the fields from the paths to the new field info + !! + subroutine updateFields(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: i + type(dictionary) :: dict + character(nameLen) :: fieldName + class(field), pointer :: ptr + class(pieceConstantField), pointer :: piecePtr + real(defReal) :: maxTemp + character(100), parameter :: Here = 'updateFields (couplingAdmin_class.f90)' + + do i = 1, size(self % fieldPaths) + call fileToDict(dict, self % fieldPaths(i)) + + ! Obtain the pointer to the corresponding field. + ! Reinitialise the field from the dictionary. + select case(trim(self % fieldNames(i))) + case('temperature') + fieldName = nameTemperature + + case('density') + fieldName = nameDensity + + case default + call fatalError(Here, 'Unrecognised field type requested') + end select + + if (.not. gr_hasField(fieldName)) then + call new_field(dict, fieldName) + end if + + ptr => gr_fieldPtr(gr_fieldIdx(fieldName)) + + call ptr % kill() + call ptr % init(dict) + + ! Ensure assumed maximum temperature is not violated + ! For now: down cast to pieceConstantField + ! TODO: implement getMaxValue for all fields + if (fieldName == nameTemperature) then + piecePtr => pieceConstantField_CptrCast(ptr) + maxTemp = piecePtr % getMaxValue() + call self % checkMaxTemperature(maxTemp) + end if + + end do + + end subroutine updateFields + + !! + !! Writes a file to signal that the current iteration is finished + !! + subroutine signalIterationOver(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: unit + + open(newunit=unit, file=self % commFileOut, status="replace", action="write") + write(unit, '(A)') CONTINUE_SIGNAL + close(unit) + + end subroutine signalIterationOver + + !! + !! Writes a file to signal that the entire calculation is finished + !! + subroutine signalCalculationOver(self) + class(couplingAdmin), intent(in) :: self + integer(shortInt) :: unit + + open(newunit=unit, file=self % commFileOut, status="replace", action="write") + write(unit, '(A)') END_SIGNAL + close(unit) + + end subroutine signalCalculationOver + + !! + !! Pauses the calculation while waiting for the existence + !! of a signal file. + !! + subroutine waitForSignal(self) + class(couplingAdmin), intent(inout) :: self + logical(defBool) :: exists + integer(shortInt) :: unit, ios + character(nameLen) :: line + character(100), parameter :: Here ='waitForSignal (couplingAdmin_class.f90)' + + do + + ! Check for the presence of a signal file + inquire(file=self % commFileIn, exist=exists) + + if (exists) then + ! Read and delete the communication file + open(newunit=unit, file=self % commFileIn, status="old") + + read(unit, '(A)', iostat=ios) line + + ! Make sure the file isn't partially written: + if (ios /= 0) then + call c_usleep(int(0.005 * 1E6)) + cycle + end if + + ! Check whether to continue iteration or end coupling + if (trim(adjustl(line)) == END_SIGNAL) then + self % isCoupled = .false. + elseif (trim(adjustl(line)) /= CONTINUE_SIGNAL) then + call fatalError(Here, 'Unrecognised signal: '//trim(adjustl(line))) + end if + + close(unit, status="delete") + return + end if + + call c_usleep(int(0.005 * 1E6)) + + end do + + end subroutine waitForSignal + + !! + !! Attach coupling tallies to another tally + !! + subroutine attachTally(self, tallyPtr) + class(couplingAdmin), intent(in) :: self + type(tallyAdmin), pointer, intent(inout) :: tallyPtr + + if (self % doCoupling()) then + call tallyPtr % push(self % tally) + end if + + end subroutine attachTally + + !! + !! Output tally results + !! + subroutine outputTallies(self) + class(couplingAdmin), intent(in) :: self + type(outputFile) :: out + + call out % init(self % outputFormat, filename = self % outputFile) + call self % tally % print(out) + + end subroutine outputTallies + +end module couplingAdmin_class diff --git a/DataStructures/charMap_class.f90 b/DataStructures/charMap_class.f90 index b1e14f798..a9b5103fe 100644 --- a/DataStructures/charMap_class.f90 +++ b/DataStructures/charMap_class.f90 @@ -34,7 +34,7 @@ module charMap_class !! Implementation is based on intMap !! !! NOTE: Following structure can be used to loop over entire map - !! But note that you CANNOT modifay entries with "add" inside this loop. + !! But note that you CANNOT modify entries with "add" inside this loop. !! Use "atSet" instead !! it = map % begin() !! do while (it /= map % end()) diff --git a/DataStructures/heapQueue_class.f90 b/DataStructures/heapQueue_class.f90 index 6daf01854..cbdcdf818 100644 --- a/DataStructures/heapQueue_class.f90 +++ b/DataStructures/heapQueue_class.f90 @@ -138,8 +138,10 @@ subroutine replace(self, val) ! We need to consider the case where there is only one child in a node ! (can happen when the size is even). If child is the last element it is the ! node with no sibling - if (child /= self % size .and. self % heap(child) < self % heap(child + 1)) then - child = child + 1 + if (child /= self % size) then + if (self % heap(child) < self % heap(child + 1)) then + child = child + 1 + end if end if ! If the parent is larger than the larger child we are done diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 index e1b109513..fd4335f1f 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 @@ -130,7 +130,7 @@ subroutine init_dict(self, dict) ! Check for invalid pitch if (any(self % pitch < 10 * SURF_TOL)) then call fatalError(Here, 'Pitch size must be larger than: '//numToChar( 10 * SURF_TOL)) - end if + end if ! Calculate halfwidth and corner self % a_bar = self % pitch * HALF - SURF_TOL @@ -163,15 +163,21 @@ subroutine init_dict(self, dict) end if ! Size field value array - self % N = product(self % sizeN * self % nMat) + self % N = product(self % sizeN) * self % nMat allocate(self % val(self % N + 1)) ! Read field values for each material idx0 = 0 do i = 1, size(mats) + if (allocated(temp)) deallocate(temp) call dict % get(temp, mats(i)) + ! Ensure size is as expected + if (size(temp) /= product(self % sizeN)) call fatalError(Here,& + 'Invalid input size to field. Should be '//numToChar(self % sizeN)//'. '//& + 'Was: '//numToChar(size(temp))) + ! Flip array up-down for more natural input ! Reshape into rank 2 array tempMap = reshape(temp, [self % sizeN(1), self % sizeN(2) * self % sizeN(3)]) @@ -347,6 +353,7 @@ elemental subroutine kill(self) self % nMat = 0 self % corner = ZERO self % a_bar = ZERO + if (allocated(self % matIdxs)) deallocate(self % matIdxs) call self % outline % kill() end subroutine kill diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 index 8d21708e5..3889b05fd 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 @@ -117,7 +117,7 @@ subroutine init(self, dict) class(pieceConstantField), intent(inout) :: self class(dictionary), intent(in) :: dict character(pathLen) :: path - class(dictionary), pointer :: tempDict + type(dictionary) :: tempDict if (dict % isPresent('path')) then call dict % get(path, 'path') diff --git a/Geometry/Tests/geometryStd_iTest.f90 b/Geometry/Tests/geometryStd_iTest.f90 index 028a53342..50cb535c7 100644 --- a/Geometry/Tests/geometryStd_iTest.f90 +++ b/Geometry/Tests/geometryStd_iTest.f90 @@ -8,6 +8,8 @@ module geometryStd_iTest use coord_class, only : coordList use geometryStd_class, only : geometryStd use materialMenu_mod, only : mm_init => init + use geometryReg_mod, only : gr_kill => kill + use fieldFactory_func, only : new_field use funit implicit none @@ -213,8 +215,62 @@ subroutine test_lattice_geom() @assertEqual(idx, coords % matIdx) @assertEqual(0.08_defReal, maxDist, TOL) + ! Add superimposed density field to complicate movement + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + + ! Place a particle and move it, ensuring the event is as expected + r = [-0.64_defReal, -1.18_defReal, 0.0_defReal] + u = [ONE, ZERO, ZERO] + call coords % init(r, u) + call geom % placeCoord(coords) + maxDist = 10.0_defReal + + ! Cross the field + r_ref = [-0.63_defReal, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(FIELD_EV, event) + @assertEqual(0.01_defReal, maxDist, TOL) + + ! Cross a cell rather than hit the field + maxDist = 1.0_defReal + r_ref = [ZERO, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(CROSS_EV, event) + @assertEqual(0.63_defReal, maxDist, TOL) + + ! Collision close to the field crossing + maxDist = 0.629_defReal + r_ref = [0.629_defReal, -1.18_defReal, ZERO] + u_ref = [ONE, ZERO, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(COLL_EV, event) + @assertEqual(0.629_defReal, maxDist, TOL) + + ! Boundary crossing + u = [ZERO, -ONE, ZERO] + call coords % assignDirection(u) + maxDist = 2.0_defReal + + r_ref = [0.629_defReal, 1.26_defReal, ZERO] + u_ref = [ZERO, -ONE, ZERO] + call geom % move(coords, maxDist, event) + @assertEqual(r_ref, coords % lvl(1) % r, TOL) + @assertEqual(u_ref, coords % lvl(1) % dir, TOL) + @assertEqual(BOUNDARY_EV, event) + @assertEqual(0.08_defReal, maxDist, TOL) + ! Kill geometry call geom % kill() + call gr_kill() end subroutine test_lattice_geom diff --git a/Geometry/geometryReg_mod.f90 b/Geometry/geometryReg_mod.f90 index a4ceb805b..f9107cdc8 100644 --- a/Geometry/geometryReg_mod.f90 +++ b/Geometry/geometryReg_mod.f90 @@ -43,6 +43,8 @@ module geometryReg_mod implicit none private + + integer(shortInt), parameter :: NO_GEOM = -7, NO_FIELD = -8 !! !! Small type to hold a single polymorphic geometry @@ -86,8 +88,8 @@ module geometryReg_mod type(charMap) :: fieldNameMap ! Save indices for speed. Avoids repeated calls to the character map - integer(shortInt) :: temperatureIdx = -8 - integer(shortInt) :: densityIdx = -8 + integer(shortInt) :: temperatureIdx = NO_FIELD + integer(shortInt) :: densityIdx = NO_FIELD contains @@ -107,13 +109,12 @@ subroutine addGeom(geom, name) class(geometry), allocatable, intent(inout) :: geom character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -7 character(100), parameter :: Here = 'addGeom (geometryReg_mod.f90)' ! Get free index - idx = geometryNameMap % getOrDefault(name, NOT_PRESENT) + idx = geometryNameMap % getOrDefault(name, NO_GEOM) - if (idx /= NOT_PRESENT) then + if (idx /= NO_GEOM) then call fatalError(Here, 'Geometry with name: '//trim(name)//' has already been defined with & &idx: '//numToChar(idx)) else @@ -145,11 +146,10 @@ end subroutine addGeom function geomIdx(name) result(idx) character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'geomIdx (geometryReg_mod.f90)' - idx = geometryNameMap % getOrDefault(name, NOT_PRESENT) - if (idx == NOT_PRESENT) then + idx = geometryNameMap % getOrDefault(name, NO_GEOM) + if (idx == NO_GEOM) then call fatalError(Here, 'Geometry: '//trim(name)//' was not found.') end if @@ -217,13 +217,12 @@ subroutine addField(kentta, name) class(field), allocatable, intent(inout) :: kentta character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -7 character(100), parameter :: Here = 'addField (geometryReg_mod.f90)' ! Get free index - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) + idx = fieldNameMap % getOrDefault(name, NO_FIELD) - if (idx /= NOT_PRESENT) then + if (idx /= NO_FIELD) then call fatalError(Here, 'Field with name: '//trim(name)//' has already been defined with & &idx: '//numToChar(idx)) else @@ -241,7 +240,7 @@ subroutine addField(kentta, name) elseif (name == nameDensity) then densityIdx = idx end if - + ! Point field call move_alloc(kentta, fields(idx) % kentta) @@ -259,7 +258,6 @@ end subroutine addField pure function hasField(name) result(exists) character(nameLen), intent(in) :: name logical(defBool) :: exists - integer(shortInt), parameter :: NOT_PRESENT = -8 integer(shortInt) :: idx select case(name) @@ -268,10 +266,10 @@ pure function hasField(name) result(exists) case(nameDensity) idx = densityIdx case default - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) + idx = fieldNameMap % getOrDefault(name, NO_FIELD) end select - exists = (idx /= NOT_PRESENT) + exists = (idx /= NO_FIELD) end function hasField @@ -290,11 +288,10 @@ end function hasField function fieldIdx(name) result(idx) character(nameLen), intent(in) :: name integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'fieldIdx (geometryReg_mod.f90)' - idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) - if (idx == NOT_PRESENT) then + idx = fieldNameMap % getOrDefault(name, NO_FIELD) + if (idx == NO_FIELD) then call fatalError(Here, 'Field: '//trim(name)//' was not found.') end if @@ -319,7 +316,6 @@ function fieldPtrName(name) result(ptr) character(nameLen), intent(in) :: name class(field), pointer :: ptr integer(shortInt) :: idx - integer(shortInt), parameter :: NOT_PRESENT = -8 character(*), parameter :: Here = 'fieldPtrName (geometryReg_mod.f90)' select case(name) @@ -328,7 +324,7 @@ function fieldPtrName(name) result(ptr) case(nameDensity) idx = densityIdx case default - idx = NOT_PRESENT + idx = NO_FIELD end select if (idx < 1 .or. idx > fieldTop) then @@ -394,6 +390,9 @@ subroutine kill() deallocate(fields) end if fieldTop = 0 + + temperatureIdx = NO_FIELD + densityIdx = NO_FIELD end subroutine kill diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index 46aec37bf..b3feaeaea 100644 --- a/Geometry/geometryStd_class.f90 +++ b/Geometry/geometryStd_class.f90 @@ -239,7 +239,7 @@ subroutine move_noCache(self, coords, maxDist, event) event = COLL_EV maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway - else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + else if (fieldDist < dist - NUDGE) then ! Stays within the same cell, but crosses field boundary call coords % moveLocal(fieldDist, level0) event = FIELD_EV maxDist = fieldDist @@ -312,7 +312,7 @@ subroutine move_withCache(self, coords, maxDist, event, cache) maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway cache % lvl = 0 - else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + else if (fieldDist < dist - NUDGE) then ! Stays within the same cell, but crosses field boundary call coords % moveLocal(fieldDist, level0) event = FIELD_EV maxDist = fieldDist diff --git a/Geometry/geometry_inter.f90 b/Geometry/geometry_inter.f90 index c473f5194..8e9161e14 100644 --- a/Geometry/geometry_inter.f90 +++ b/Geometry/geometry_inter.f90 @@ -1,7 +1,7 @@ module geometry_inter use numPrecision - use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INFINITY, COLL_EV + use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INF, COLL_EV use genericProcedures, only : fatalError use dictionary_class, only : dictionary use charMap_class, only : charMap @@ -552,14 +552,14 @@ subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) ! Trace until an opaque material is hit do while (any(mats == matIdx)) - dist = INFINITY + dist = INF call self % moveNoBC(ray, dist, event, normal0) matIdx = ray % matIdx ! If distance is infinite or particle collided, not pointing towards anything opaque. ! Hence, short circuit the trace - if ((dist == INFINITY) .or. (event == COLL_EV)) then + if ((dist == INF) .or. (event == COLL_EV)) then lightTrace = .false. exit end if diff --git a/InputFiles/Benchmarks/Transients/LatticeBrunner b/InputFiles/Benchmarks/Transients/LatticeBrunner index bcdb236ab..c357d10db 100644 --- a/InputFiles/Benchmarks/Transients/LatticeBrunner +++ b/InputFiles/Benchmarks/Transients/LatticeBrunner @@ -8,12 +8,12 @@ type kineticPhysicsPackage; -pop 5000000; +pop 500000; timeSteps 1; XSdata mg; dataType mg; dt 7.0; -cycles 750; +cycles 250; buffer 100; collisionOperator { neutronMG {type neutronMGstd;}} @@ -109,7 +109,7 @@ nuclearData { temp 12345; composition { } - xsFile ./XS/abs; + xsFile ./XS/air; } diff --git a/IntegrationTestFiles/Coupling/test_coupling b/IntegrationTestFiles/Coupling/test_coupling new file mode 100644 index 000000000..703382c07 --- /dev/null +++ b/IntegrationTestFiles/Coupling/test_coupling @@ -0,0 +1,19 @@ +tally { + norm totalPow; + normVal 70000; + totalPow { type collisionClerk; response (fission); fission { type macroResponse; MT -80;}} + axialPow {type collisionClerk; + map {type spaceMap; axis z; grid lin; min -200; max 200.0; N 40;} + response (fission); fission { type macroResponse; MT -80;} + } +} + +maxTemp 1100.0; +receiveFile ./IntegrationTestFiles/Coupling/toSCONE.dat; +sendFile ./IntegrationTestFiles/Coupling/toDriver.dat; +outputFile ./IntegrationTestFiles/Coupling/sconeOut; +updateFreq 5; +maxIt 10; +outputFormat asciiMATLAB; +fieldPaths (./IntegrationTestFiles/temp ./IntegrationTestFiles/dens); +fieldNames (temperature density); diff --git a/IntegrationTestFiles/Geometry/test_lat b/IntegrationTestFiles/Geometry/test_lat index df07f69ff..96a93bd9e 100644 --- a/IntegrationTestFiles/Geometry/test_lat +++ b/IntegrationTestFiles/Geometry/test_lat @@ -45,5 +45,12 @@ nuclearData { } } - - +density { +type cartesianField; +materials (all); +origin (0 0 0); +shape (4 4 0); +pitch (0.63 0.63 0); +default -1; +all ( 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1); +} diff --git a/IntegrationTestFiles/dens b/IntegrationTestFiles/dens new file mode 100644 index 000000000..8670240f7 --- /dev/null +++ b/IntegrationTestFiles/dens @@ -0,0 +1,7 @@ +type cartesianField; +materials (cool); +origin (0 0 0); +shape (1 1 0); +pitch (1.26 1.26 0); +default -1; +cool ( 0.9 ); diff --git a/IntegrationTestFiles/temp b/IntegrationTestFiles/temp new file mode 100644 index 000000000..d53eaac1d --- /dev/null +++ b/IntegrationTestFiles/temp @@ -0,0 +1,7 @@ +type cartesianField; +materials (fuel); +origin (0 0 0); +shape (1 1 0); +pitch (1.26 1.26 0); +default -1; +fuel ( 900 ); diff --git a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 index e569690bc..7829910a4 100644 --- a/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 +++ b/NuclearData/Reactions/thermalScattReactionCE/thermalScatterInelastic_class.f90 @@ -2,6 +2,7 @@ module thermalScatterInelastic_class use numPrecision use endfConstants + use universalVariables, only : MAXIMUM_ENERGY, MINIMUM_ENERGY use genericProcedures, only : binarySearch, fatalError, numToChar, endfInterpolate, searchError use RNG_class, only : RNG use dataDeck_inter, only : dataDeck @@ -295,10 +296,13 @@ subroutine sampleOut(self, mu, phi, E_out, E_in, rand, lambda) end if - if (E_out > 20.0 .or. E_out <= ZERO) then + if (E_out > MAXIMUM_ENERGY .or. E_out <= ZERO) then call fatalError(Here,'Failed to find energy: '//numToChar(E_out)) end if + ! Allow energies lower than 1E-11 to be sampled, but floor to the minimum + E_out = max(E_out, MINIMUM_ENERGY) + ! Sample phi phi = rand % get() * TWO_PI diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 index a16ec7a1c..5143dc950 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/thermalScatteringData_iTest.f90 @@ -80,9 +80,9 @@ subroutine test_thermalScatteringData() ! Get nuclides O16 => aceNeutronNuclide_CptrCast( data % getNuclide(1)) - H1 => aceNeutronNuclide_CptrCast( data % getNuclide(2)) + C12 => aceNeutronNuclide_CptrCast( data % getNuclide(2)) H1_2 => aceNeutronNuclide_CptrCast( data % getNuclide(3)) - C12 => aceNeutronNuclide_CptrCast( data % getNuclide(4)) + H1 => aceNeutronNuclide_CptrCast( data % getNuclide(4)) !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test scattering tables @@ -158,8 +158,8 @@ subroutine test_thermalScatteringData() !<><><><><><><><><><><><><><><><><><><><><><><><> ! Test getting XSs ! H-1 - nuc => ceNeutronNuclide_CptrCast(data % getNuclide(2)) - nuclideCache(2) % E_tot = ONE + nuc => ceNeutronNuclide_CptrCast(data % getNuclide(4)) + nuclideCache(4) % E_tot = ONE call nuc % getMicroXSs(microXSs, 1.8E-6_defReal, ZERO, p % pRNG) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 8492876d5..377a6ebf2 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -10,6 +10,7 @@ module aceNeutronDatabase_class use RNG_class, only : RNG use charMap_class, only : charMap use intMap_class, only : intMap + use hashFunctions_func, only : FNV_1 ! Nuclear Data Interfaces use nuclearDatabase_inter, only : nuclearDatabase @@ -123,6 +124,7 @@ module aceNeutronDatabase_class procedure :: updateTotalTempMajXS procedure :: updateRelEnMacroXSs procedure :: makeNuclideName + procedure :: computeMaxTrackXS end type aceNeutronDatabase @@ -416,14 +418,14 @@ subroutine updateTrackMatXS(self, E, matIdx, temp, rho, rand) if (present(temp)) then T = temp else - T = -ONE + T = NO_TEMPERATURE end if matCache % T_track = T if (present(rho)) then densityFactor = rho else - densityFactor = -ONE + densityFactor = NO_DENSITY end if matCache % rho_track = densityFactor @@ -452,6 +454,8 @@ end subroutine updateTrackMatXS !! Args: !! E [in] -> Incident neutron energy for which temperature majorant is found !! matIdx [in] -> Index of material for which the material temperature majorant is found + !! temp [in] -> Temperature, in Kelvin, locally + !! rho [in] -> Density scaling factor (dimensionless) locally !! subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) class(aceNeutronDatabase), intent(in) :: self @@ -469,7 +473,7 @@ subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) matCache % trackXS = ZERO ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -499,6 +503,52 @@ subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) end associate end subroutine updateTotalTempMajXS + + !! + !! This subroutine exists to simplify the construction of the majorant in the presence + !! of a temperature field. It updates the materialCache to have the largest cross section + !! of two if evaluated at two different temperatures: an imposed temperature and the input + !! material temperature. + !! + !! The majorant cross section at a given energy point may be given by either the higher + !! or lower of the possible temperatures. If a temperature field is imposed (where the + !! temperature will vary spatially across a given material) then the majorant should + !! be chosen appropriately/conservatively. + !! + subroutine computeMaxTrackXS(self, E, matIdx, maxTemp, rho) + class(aceNeutronDatabase), intent(in) :: self + real(defReal), intent(in) :: E + integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: maxTemp + real(defReal), intent(in) :: rho + real(defReal) :: xs1, xs2 + + associate (matCache => cache_materialCache(matIdx)) + + ! Clean current total XS + matCache % trackXS = ZERO + + ! If no temperature is given then the usual majorant from the material + ! temperature is fine and we needn't go further + if (maxTemp < ZERO) then + call self % updateTotalTempMajXS(E, matIdx, maxTemp, rho) + + ! Otherwise, compute the xs for both the material temperature and the given + ! max temperature + else + call self % updateTotalTempMajXS(E, matIdx, NO_TEMPERATURE, rho) + xs1 = matCache % trackXS + + call self % updateTotalTempMajXS(E, matIdx, maxTemp, rho) + xs2 = matCache % trackXS + + matCache % trackXS = max(xs1, xs2) + + end if + + end associate + + end subroutine computeMaxTrackXS !! !! Make sure that totalXS of material with matIdx is at energy E @@ -525,7 +575,7 @@ subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) matCache % rho_tot = rho ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -563,7 +613,6 @@ subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) matCache % xss % total = matCache % xss % total * densityFactor end if - matCache % T_tot = temp end associate @@ -594,7 +643,7 @@ subroutine updateMacroXSs(self, E, matIdx, temp, rho, rand) matCache % rho_tot = rho ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = mat % kT else kT = temp * kBoltzmannMeV @@ -832,7 +881,7 @@ subroutine updateTotalTempNucXS(self, E, kT, nucIdx) ! Catch negative dkT - perhaps due to input field values if (deltakT < ZERO) call fatalError(Here,& - 'Invalid input temperatureL '//numToChar(kT/kBoltzmannMeV)) + 'Invalid input temperature: '//numToChar(kT/kBoltzmannMeV)) ! Check if an update is required if (nucCache % E_maj /= E .or. nucCache % deltakT /= deltakT) then @@ -879,19 +928,21 @@ subroutine init(self, dict, ptr, silent ) type(materialItem), pointer :: mat class(ceNeutronDatabase), pointer :: ptr_ceDatabase type(charMap) :: nucSet + type(dictionary) :: sabDict type(aceCard) :: ACE - type(aceSabCard) :: ACE_Sab1, ACE_Sab2 + type(aceSabCard), dimension(:), allocatable :: ACE_Sab + character(nameLen), dimension(:), allocatable :: fileNames character(pathLen) :: aceLibPath - character(nameLen) :: name, name_file1, name_file2, nucDBRC_temp - integer(shortInt) :: i, j, envFlag, nucIdx, idx, idx1, idx2 + character(nameLen) :: nucDBRC_temp, name + integer(shortInt) :: i, j, envFlag, nucIdx, idx integer(shortInt) :: maxNuc logical(defBool) :: isFissileMat integer(shortInt),dimension(:),allocatable :: nucIdxs, zaidDBRC character(nameLen),dimension(:),allocatable :: nucDBRC real(defReal) :: A, nuckT, eUpSab, eUpSabNuc, & eLowURR, eLowUrrNuc, alpha, & - deltakT, eUpper, eLower, kT, & - temp + deltakT, eUpper, eLower, & + temp, kT real(defReal), dimension(2) :: sabT integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 character(100), parameter :: Here = 'init (aceNeutronDatabase_class.f90)' @@ -912,6 +963,9 @@ subroutine init(self, dict, ptr, silent ) ptr_ceDatabase => ceNeutronDatabase_CptrCast(ptr) if(.not.associated(ptr_ceDatabase)) call fatalError(Here,"Should not happen. WTF?!") + ! Initialise dictionary to store S(a,b) info + call sabDict % init(1) + ! Create list of all nuclides. Loop over materials ! Find maximum number of nuclides: maxNuc do i = 1, mm_nMat() @@ -922,6 +976,11 @@ subroutine init(self, dict, ptr, silent ) do j = 1, size(mat % nuclides) name = self % makeNuclideName(mat % nuclides(j)) call nucSet % add(name, IN_SET) + + ! Add S(alpha,beta) data + if (mat % nuclides(j) % hasSab) then + call sabDict % store(name, mat % nuclides(j) % file_Sab) + end if end do end do @@ -994,41 +1053,41 @@ subroutine init(self, dict, ptr, silent ) nucIdx = 1 do while (i /= nucSet % end()) - idx1 = index(nucSet % atKey(i),'+') - idx2 = index(nucSet % atKey(i),'#') - if (idx1 /= 0) then - name = trim(nucSet % atKey(i)) - if (idx2 == 0) then - name_file1 = trim(name(idx1+1:nameLen)) - else - name_file1 = trim(name(idx1+1:idx2-1)) - name_file2 = trim(name(idx2+1:nameLen)) - end if - name = name(1:idx1-1) - else - name = nucSet % atKey(i) + name = trim(nucSet % atKey(i)) + + ! Get Sab filenames + if (sabDict % isPresent(name)) then + call sabDict % get(fileNames, name) + + ! Truncate name to remove hash + idx = index(name,'_') + if (idx == 0) call fatalError(Here,'Failed to parse Sab nuclide: '//name) + name = name(1:idx-1) end if if (loud) then call statusMsg("Building: "// trim(name)// " with index: " //numToChar(nucIdx)) - if (idx1 /= 0 .and. idx2 == 0) & - call statusMsg("including S(alpha,beta) tables with file: " //trim(name_file1)) - if (idx1 /= 0 .and. idx2 /= 0) & - call statusMsg("including S(alpha,beta) tables with files: " //trim(name_file1)//' '//trim(name_file2)) + if (allocated(fileNames)) then + call statusMsg("including S(alpha,beta) tables with files: ") + do idx = 1, size(fileNames) + call statusMsg(fileNames(idx)) + end do + end if end if call new_neutronACE(ACE, name) call self % nuclides(nucIdx) % init(ACE, nucIdx, ptr_ceDatabase) ! Initialise S(alpha,beta) tables - if (idx1 /= 0 ) then - call new_moderACE(ACE_Sab1, name_file1) - if (idx2 /= 0) then - call new_moderACE(ACE_Sab2, name_file2) - call self % nuclides(nucIdx) % initSab(ACE_Sab1, ACE_Sab2) - else - call self % nuclides(nucIdx) % initSab(ACE_Sab1) - end if + if (allocated(fileNames)) then + allocate(ACE_Sab(size(fileNames))) + do idx = 1, size(fileNames) + call new_moderACE(ACE_Sab(idx), fileNames(idx)) + end do + call self % nuclides(nucIdx) % initSab(ACE_Sab) + fileNames = '' + deallocate(fileNames) + deallocate(ACE_Sab) end if ! Initialise probability tables @@ -1077,7 +1136,7 @@ subroutine init(self, dict, ptr, silent ) kT = mat % T * kBoltzmannMev if ((kT < sabT(1)) .or. (kT > sabT(2))) call fatalError(Here,& 'Material temperature must be bounded by the provided S(alpha,beta) data. '//& - 'The material temperature is '//numToChar(kT / kBoltzmannMeV)//& + 'The material temperature is '//numToChar(mat % T)//& 'K while the data bounds are '//numToChar(sabT(1) / kBoltzmannMeV)//& 'K and '//numToChar(sabT(2) / kBoltzmannMeV)//'K.') end if @@ -1171,7 +1230,8 @@ function makeNuclideName(self, nuclide) result(name) class(aceNeutronDatabase), intent(in) :: self type(nuclideInfo), intent(in) :: nuclide character(nameLen) :: name - character(:), allocatable :: file + character(:), allocatable :: extension + integer(shortInt) :: i, hash name = trim(nuclide % toChar()) @@ -1180,17 +1240,13 @@ function makeNuclideName(self, nuclide) result(name) ! scattering if (nuclide % hasSab) then - file = trim(nuclide % file_Sab1) - name = trim(name) // '+' // file - deallocate(file) - - ! Attach second Sab file for stochastic mixing - if (nuclide % sabMix) then - file = trim(nuclide % file_Sab2) - name = trim(name) // '#' // file - deallocate(file) - end if + extension = '' + do i = 1, size(nuclide % file_Sab) + extension = trim(extension)//'_'// trim(nuclide % file_Sab(i)) + end do + call FNV_1(extension, hash) + name = trim(name) //'_'// numToChar(hash) end if end function makeNuclideName @@ -1353,12 +1409,12 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) ! Check maxTemp is present and sensible if (present(maxTemp)) then if (maxTemp < ZERO) then - broadenTemp = -ONE + broadenTemp = NO_TEMPERATURE else broadenTemp = maxTemp end if else - broadenTemp = -ONE + broadenTemp = NO_TEMPERATURE end if ! Check scaleDensity is present and sensible @@ -1568,8 +1624,9 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) matIdx = self % activeMat(j) ! Get material tracking cross section - if (self % materials(matIdx) % hasTMS .and. broadenTemp > ZERO) then - call self % updateTrackMatXS(E, matIdx, temp = broadenTemp, rho = densityFactor, rand = rand) + if (self % materials(matIdx) % useTMS(E) .and. broadenTemp > ZERO) then + call self % computeMaxTrackXS(E, matIdx, maxTemp = broadenTemp, rho = densityFactor) + !call self % updateTrackMatXS(E, matIdx, temp = broadenTemp, rho = densityFactor, rand = rand) else call self % updateTrackMatXS(E, matIdx, rho = densityFactor, rand = rand) end if @@ -1597,8 +1654,8 @@ subroutine initMajorant(self, loud, maxTemp, scaleDensity) nucXS = urrMaj end if - ! Update total material cross section - trackXS = trackXS + dens * (nucXS - cache_nuclideCache(nucIdx) % xss % total) + ! Update total material cross section + trackXS = trackXS + dens * (nucXS - cache_nuclideCache(nucIdx) % xss % total) end if diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 index 2a7a5377f..23654f2c3 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronNuclide_class.f90 @@ -322,6 +322,24 @@ elemental subroutine kill(self) if (allocated(self % mainData)) deallocate(self % mainData) if (allocated(self % eGrid)) deallocate(self % eGrid) call self % idxMT % kill() + + ! Clean S(alpha, beta) + if (allocated(self % thData)) then + do i = 1, size(self % thData) + call self % thData(i) % kill() + end do + deallocate(self % thData) + self % SabEl = ZERO + self % SabInel = ZERO + self % hasThData = .false. + self % stochasticMixing = .false. + end if + + ! Clean URR + self % urrE = ZERO + call self % probTab % kill() + self % hasProbTab = .false. + self % IFF = 0 end subroutine kill @@ -988,67 +1006,85 @@ subroutine initUrr(self, ACE) end subroutine initUrr !! - !! Initialise thermal scattering tables from ACE card + !! Initialise thermal scattering tables from ACE card. + !! If the array has size greater than 1, activates stochastic mixing treatment. !! !! Args: - !! ACE1 [inout] -> ACE S(a,b) card - !! ACE2 [inout] -> Optional second ACE S(a,b) card + !! ACE [inout] -> array of ACE S(a,b) cards !! !! Errors: !! fatalError if the inelastic scattering S(a,b) energy grid starts at a - !! lower energy than the nuclide energy grid - !! - subroutine initSab(self, ACE1, ACE2) - class(aceNeutronNuclide), intent(inout) :: self - class(aceSabCard), intent(inout) :: ACE1 - class(aceSabCard), intent(inout), optional :: ACE2 - real(defReal), dimension(2) :: EBounds - real(defReal) :: T1, T2 - type(thermalData) :: temp + !! lower energy than the nuclide energy grid. + !! + !! fatalError if any of the input temperatures are equal. + !! + subroutine initSab(self, ACE) + class(aceNeutronNuclide), intent(inout) :: self + class(aceSabCard), dimension(:), intent(inout) :: ACE + real(defReal), dimension(2) :: EBounds + real(defReal) :: T + integer(shortInt) :: nFiles + type(thermalData) :: temp + integer(shortInt) :: i, j character(100), parameter :: Here = "initSab (aceNeutronNuclide_class.f90)" - if (present(ACE2)) then - allocate(self % thData(2)) - else - allocate(self % thData(1)) - end if - - ! Initialise S(a,b) class from ACE file - call self % thData(1) % init(ACE1) self % hasThData = .true. + nFiles = size(ACE) + allocate(self % thData(nFiles)) + if (nFiles > 1) self % stochasticMixing = .true. - ! Initialise energy boundaries - self % SabInel = self % thData(1) % getEBounds('inelastic') - self % SabEl = self % thData(1) % getEBounds('elastic') + ! Initialise S(a,b) class from ACE file + do i = 1, nFiles + call self % thData(i) % init(ACE(i)) + + if (i == 1) then + ! Initialise energy boundaries + self % SabInel = self % thData(1) % getEBounds('inelastic') + self % SabEl = self % thData(1) % getEBounds('elastic') - ! Add second S(a,b) file for stochastic mixing - if (present(ACE2)) then + else + ! Ensure energy bounds are conservative + EBounds = self % thData(i) % getEBounds('inelastic') + if (EBounds(1) > self % SabInel(1)) self % SabInel(1) = EBounds(1) + if (EBounds(2) < self % SabInel(2)) self % SabInel(2) = EBounds(2) + + EBounds = self % thData(i) % getEbounds('elastic') + if (EBounds(1) > self % SabEl(1)) self % SabEl(1) = EBounds(1) + if (EBounds(2) < self % SabEl(2)) self % SabEl(2) = EBounds(2) - self % stochasticMixing = .true. - call self % thData(2) % init(ACE2) + end if - ! Ensure energy bounds are conservative - EBounds = self % thData(2) % getEBounds('inelastic') - if (EBounds(1) > self % SabInel(1)) self % SabInel(1) = EBounds(1) - if (EBounds(2) < self % SabInel(2)) self % SabInel(2) = EBounds(2) + end do - EBounds = self % thData(2) % getEbounds('elastic') - if (EBounds(1) > self % SabEl(1)) self % SabEl(1) = EBounds(1) - if (EBounds(2) < self % SabEl(2)) self % SabEl(2) = EBounds(2) + ! Rearrange files to be in order of ascending temperature + if (nFiles > 1) then + do i = 2, nFiles + temp = self % thData(i) + T = temp % getTemperature() + j = i - 1 - ! Identify which data is higher temperature and which is lower - ! 1 should be lower than 2 - swap if necessary - T1 = self % thData(1) % getTemperature() - T2 = self % thData(2) % getTemperature() + do while (j >= 1) + if (self % thData(j) % getTemperature() <= T) exit + self % thData(j+1) = self % thData(j) + j = j - 1 + end do - if (T1 > T2) then - temp = self % thData(1) - self % thData(1) = self % thData(2) - self % thData(2) = temp - end if + self % thData(j+1) = temp + end do + + ! Check to ensure all temperatures are unique + T = self % thData(1) % getTemperature() / kBoltzmannMeV + do i = 2, nFiles + if (abs(T - self % thData(i) % getTemperature()) / kBoltzmannMeV < 1.0E-6) call fatalError(& + Here,'S(alpha,beta) files have been input with the same temperature. '//& + 'Data '//numToChar(i-1)//': '//numToChar(T)//'. Data '//numToChar(i)//& + ': '//numToChar(self % thData(i) % getTemperature()/kBoltzmannMeV)) + T = self % thData(i) % getTemperature() / kBoltzmannMeV + end do end if + ! Check consistency of energy grid if (self % SabInel(1) < self % eGrid(1)) then call fatalError(Here, 'S(a,b) low energy boundary is lower than nuclide first energy point') @@ -1064,6 +1100,11 @@ end subroutine initSab !! Returns the Sab index for later consistent handling !! of angular distributions by storing in a cache !! + !! Errors: + !! fatalError if requested temperature is outside of the possible temperature range + !! + !! fatalError if the linear search over temperatures for interpolation fails + !! subroutine getSabPointer(self, kT, rand, ptr, idx) class(aceNeutronNuclide), intent(in), target :: self real(defReal), intent(in) :: kT @@ -1071,23 +1112,35 @@ subroutine getSabPointer(self, kT, rand, ptr, idx) type(thermalData), pointer, intent(out) :: ptr integer(shortInt), intent(out) :: idx real(defReal) :: kT1, kT2 + real(defReal), dimension(2) :: kTMax + integer(shortInt) :: i character(100), parameter :: Here = "getSabPointer (aceNeutronNuclide_class.f90)" if (self % stochasticMixing) then - kT1 = self % thData(1) % getTemperature() - kT2 = self % thData(2) % getTemperature() + kTMax = self % getSabTBounds() - if ((kT < kT1) .or. (kT > kT2)) call fatalError(Here,& + if ((kT < kTMax(1)) .or. (kT > kTMax(2))) call fatalError(Here,& 'Requested temperature '//numToChar(kT)//' not in temperature bounds: '//& - numToChar(kT1)//' and '//numToChar(kT2)) + numToChar(kTMax(1))//' and '//numToChar(kTMax(2))) + + ! Search through temperatures to find the interpolation range + do i = 1, size(self % thData) - 1 + if (kT <= self % thData(i+1) % getTemperature()) then + kT1 = self % thData(i) % getTemperature() + kT2 = self % thData(i+1) % getTemperature() + + if ((kT2 - kT)/(kT2 - kT1) > rand % get()) then + idx = i + else + idx = i + 1 + end if + ptr => self % thData(idx) + return + end if + end do - if ((kT2 - kT)/(kT2 - kT1) > rand % get()) then - ptr => self % thData(1) - idx = 1 - else - ptr => self % thData(2) - idx = 2 - end if + call fatalError(Here,'Temperature search failed. Requested temperature '//numToChar(kT)//& + ' exceed maximum: '//numToChar(self % thData(i+1) % getTemperature())) else ptr => self % thData(1) idx = 1 @@ -1099,13 +1152,13 @@ end subroutine getSabPointer !! Return the temperature bounds of S(alpha,beta) data !! If only one library, both bounds are the same temperature !! - function getSabTBounds(self) result(kT) + pure function getSabTBounds(self) result(kT) class(aceNeutronNuclide), intent(in) :: self real(defReal), dimension(2) :: kT kT(1) = self % thData(1) % getTemperature() if (self % stochasticMixing) then - kT(2) = self % thData(2) % getTemperature() + kT(2) = self % thData(size(self % thData)) % getTemperature() else kT(2) = kT(1) end if @@ -1117,7 +1170,7 @@ end function getSabTBounds !! !! Prints: !! nucIdx; Size of Energy grid; Maximum and Minumum energy; Mass; Temperature; ACE ZAID; - !! MT numbers used in tranposrt (active MTs) and not used in transport (inactiveMTs) + !! MT numbers used in transport (active MTs) and not used in transport (inactiveMTs) !! !! NOTE: For Now Formatting is horrible. Can be used only for debug !! MT = 18 may appear as inactive MT but is included from FIS block ! @@ -1151,7 +1204,6 @@ subroutine display(self) allMT = size(self % MTdata) print '(A)', "Active MTs: " // numToChar(self % MTdata(1:sMT) % MT) print '(A)', "Inactive MTs: "// numToChar(self % MTdata(sMT+1:allMT) % MT) - print '(A)', "This implementation ignores MT=5 (N,anything) !" end subroutine display diff --git a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 index 60c5402c8..b3b8e7154 100644 --- a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 +++ b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 @@ -19,6 +19,7 @@ module ceNeutronCache_mod use numPrecision + use universalVariables, only : NO_TEMPERATURE, NO_DENSITY use genericProcedures, only : fatalError, numToChar use neutronXsPackages_class, only : neutronMicroXSs, neutronMacroXSs @@ -37,16 +38,16 @@ module ceNeutronCache_mod !! idx -> Index on a nuclide grid for energy E_tot !! xss -> Cached cross-section values !! E_track -> Energy of the tracking xs + !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS !! T_track -> Temperature at which tracking XS was evaluated !! rho_track -> Density scaling at which tracking XS was evaluated - !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS !! E_rel -> Base energy for which relative energy cross sections are found (for TMS) !! xssRel -> Cached effective cross-section values at energy relative to E_rel (for TMS) !! type, public :: cacheMatDat real(defReal) :: E_tot = ZERO - real(defReal) :: T_tot = -ONE - real(defReal) :: rho_tot = -ONE + real(defReal) :: T_tot = NO_TEMPERATURE + real(defReal) :: rho_tot = NO_DENSITY real(defReal) :: E_tail = ZERO real(defReal) :: f = ZERO integer(shortInt) :: idx = 0 @@ -55,8 +56,8 @@ module ceNeutronCache_mod ! Tracking data real(defReal) :: E_track = ZERO real(defReal) :: trackXS = ZERO - real(defReal) :: T_track = -ONE - real(defReal) :: rho_track = -ONE + real(defReal) :: T_track = NO_TEMPERATURE + real(defReal) :: rho_track = NO_DENSITY ! TMS data real(defReal) :: E_rel = ZERO diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index 96d2c41fd..20ef6f473 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -86,8 +86,8 @@ module ceNeutronDatabase_inter !! !! Return energy bounds for data in the database !! - !! eMin and eMax are minimun and maximumum energy such that data - !! for ALL nuclides if avalible + !! eMin and eMax are minimum and maximumum energy such that data + !! for ALL nuclides is available !! !! Args: !! eMin [out] -> minimum value of energy [MeV] @@ -477,7 +477,7 @@ function getMajorantXS(self, p) result(xs) ! Check dynamic type of the particle if (p % isMG .or. p % type /= P_NEUTRON) then - call fatalError(Here, 'Dynamic type of the partcle is not CE Neutron but:'//p % typeToChar()) + call fatalError(Here, 'Dynamic type of the particle is not CE Neutron but:'//p % typeToChar()) end if ! Check Cache and update if needed @@ -515,5 +515,4 @@ pure function ceNeutronDatabase_CptrCast(source) result(ptr) end function ceNeutronDatabase_CptrCast - end module ceNeutronDatabase_inter diff --git a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 index 05395c74f..b329eb30e 100644 --- a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 +++ b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 @@ -360,7 +360,7 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) trackMatXS = materialCache(self % matIdx) % trackXS * rand % get() ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -373,8 +373,8 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) densityFactor = rho end if - ! Loop over nuclides - do i = 1,size(self % nuclides) + ! Loop over nuclides until one is sampled + nucLoop: do i = 1,size(self % nuclides) nucIdx = self % nuclides(i) dens = self % dens(i) * densityFactor @@ -440,17 +440,17 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) end if - ! Exit function, return the sampled nucIdx return end if end associate - end do + end do nucLoop ! Print error message as the inversion failed call fatalError(Here,'Nuclide sampling loop failed to terminate') + end subroutine sampleNuclide @@ -503,7 +503,7 @@ function sampleFission(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -599,7 +599,7 @@ function sampleScatter(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV @@ -691,7 +691,7 @@ function sampleScatterWithFission(self, E, temp, rho, rand) result(nucIdx) end if ! Use imposed temperature if given - if (temp <= ZERO) then + if (temp < ZERO) then kT = self % kT else kT = temp * kBoltzmannMeV diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 3231ab5d1..4cd8b13a5 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -46,8 +46,7 @@ module materialMenu_mod !! T -> Evaluation number !! hasSab -> Does the nuclide have S(a,b) data? !! sabMix -> Does the nuclide mix S(a,b) data? - !! file_Sab1 -> First (and maybe only) S(a,b) file - !! file_Sab2 -> Second S(a,b) file + !! file_Sab -> Array of S(a,b) file names (generally >1 for stochastic interpolation) !! !! Interface: !! init -> build from a string @@ -58,8 +57,7 @@ module materialMenu_mod integer(shortInt) :: T = -1 logical(defBool) :: hasSab = .false. logical(defBool) :: sabMix = .false. - character(nameLen) :: file_Sab1 - character(nameLen) :: file_Sab2 + character(nameLen), allocatable, dimension(:) :: file_Sab contains procedure :: init => init_nuclideInfo procedure :: toChar => toChar_nuclideInfo @@ -354,17 +352,11 @@ subroutine init_materialItem(self, name, idx, dict) ! Check for stochastic mixing - this will depend on the ! size of the array of files produce call moderDict % get(filenames, keys(i)) - if (size(filenames) == 2) then - self % nuclides(i) % file_Sab1 = filenames(1) - self % nuclides(i) % file_Sab2 = filenames(2) - self % nuclides(i) % sabMix = .true. - elseif (size(filenames) == 1) then - self % nuclides(i) % file_Sab1 = filenames(1) - else - print *,filenames - call fatalError(Here,'Unexpectedly long moder contents. Should be 1 or 2 '//& - 'entries.') - end if + allocate(self % nuclides(i) % file_Sab(size(filenames))) + self % nuclides(i) % file_Sab = filenames + if (size(filenames) > 1) self % nuclides(i) % sabMix = .true. + filenames = '' + deallocate(filenames) end if ! Initialise the nuclides diff --git a/NuclearData/nuclearDataReg_mod.f90 b/NuclearData/nuclearDataReg_mod.f90 index 8b16e96e7..228204356 100644 --- a/NuclearData/nuclearDataReg_mod.f90 +++ b/NuclearData/nuclearDataReg_mod.f90 @@ -291,7 +291,7 @@ subroutine activate(type, name, activeMat, silent) if(idx == 0 ) then call fatalError(Here, trim(name)//' is was not defined. Cannot activate it!') else if(idx < 0) then - call fatalError(Here, '-ve idx from databaseNameMap. Quite immpossible. WTF?') + call fatalError(Here, '-ve idx from databaseNameMap. Quite impossible. WTF?') end if ! Make if it is not already made diff --git a/PhysicsPackages/alphaPhysicsPackage_class.f90 b/PhysicsPackages/alphaPhysicsPackage_class.f90 index 449e4d561..d21789b19 100644 --- a/PhysicsPackages/alphaPhysicsPackage_class.f90 +++ b/PhysicsPackages/alphaPhysicsPackage_class.f90 @@ -30,13 +30,15 @@ module alphaPhysicsPackage_class ! Geometry use geometry_inter, only : geometry use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx, & - gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_fieldPtrName => fieldPtrName use geometryFactory_func, only : new_geometry ! Fields use field_inter, only : field use uniFissSitesField_class, only : uniFissSitesField, uniFissSitesField_TptrCast use fieldFactory_func, only : new_field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast ! Nuclear Data use materialMenu_mod, only : mm_nMat => nMat @@ -62,6 +64,9 @@ module alphaPhysicsPackage_class use tallyResult_class, only : tallyResult use kAlphaAnalogClerk_class, only : kAlphaResult + ! Coupling + use couplingAdmin_class, only : couplingAdmin + ! Factories use transportOperatorFactory_func, only : new_transportOperator @@ -105,6 +110,9 @@ module alphaPhysicsPackage_class integer(shortInt) :: bufferSize logical(defBool) :: UFS = .false. + ! Coupling handler + type(couplingAdmin) :: couplingInfo + ! Calculation components type(particleDungeon), pointer :: thisCycle => null() type(particleDungeon), pointer :: nextCycle => null() @@ -142,6 +150,7 @@ subroutine run(self) call self % generateInitialState() call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) call self % buildActiveAtch() + call self % couplingInfo % moveToActive(self % activeTally) call self % cycles(self % activeTally, self % activeAtch, self % N_active) ! Collect results from other processes @@ -328,6 +337,9 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) call statusMsg("Time to end: " // trim(secToChar(T_toEnd))) call tally % display() + ! Perform coupling operations + call self % couplingInfo % couple(i) + end do ! Load elapsed time @@ -422,6 +434,8 @@ subroutine init(self, dict) type(outputFile) :: test_out type(visualiser) :: viz class(field), pointer :: field + class(pieceConstantField), pointer :: pcField + real(defReal) :: maxTemperature, maxDensityScale character(100), parameter :: Here ='init (alphaPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -503,6 +517,31 @@ subroutine init(self, dict) ! Activate Nuclear Data *** All materials are active call ndReg_activate(self % particleType, nucData, self % geom % activeMats()) self % nucData => ndReg_get(self % particleType) + + ! If present, build temperature field + if (dict % isPresent('temperature')) then + tempDict => dict % getDictPtr('temperature') + call new_field(tempDict, nameTemperature) + field => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(field) + maxTemperature = pcField % getMaxValue() + else + maxTemperature = NO_TEMPERATURE + end if + + ! If present, build density field + if (dict % isPresent('density')) then + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + field => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(field) + maxDensityScale = pcField % getMaxValue() + else + maxDensityScale = NO_DENSITY + end if + + ! Update majorant in case of density and temperature fields + call self % nucData % initMajorant(.false., maxTemp = maxTemperature, scaleDensity = maxDensityScale) ! Call visualisation if (dict % isPresent('viz')) then @@ -590,6 +629,11 @@ subroutine init(self, dict) ! Don't create the active tally yet as the guess value ! strongly affects the convergence rate and may bias ! the active cycle estimate + + ! Attach a tally admin for coupling + if (self % couplingInfo % doCoupling()) then + call self % couplingInfo % attachTally(self % inactiveTally) + end if call self % printSettings() diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index 76829770e..2d0f1fa5d 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -65,6 +65,9 @@ module eigenPhysicsPackage_class use tallyResult_class, only : tallyResult use keffAnalogClerk_class, only : keffResult + ! Coupling + use couplingAdmin_class, only : couplingAdmin + ! Factories use transportOperatorFactory_func, only : new_transportOperator @@ -107,6 +110,9 @@ module eigenPhysicsPackage_class integer(shortInt) :: bufferSize logical(defBool) :: UFS = .false. logical(defBool) :: reproducible = .true. + + ! Coupling handler + type(couplingAdmin) :: couplingInfo ! Calculation components type(particleDungeon), pointer :: thisCycle => null() @@ -144,8 +150,14 @@ subroutine run(self) call self % generateInitialState() call self % cycles(self % inactiveTally, self % inactiveAtch, self % N_inactive) + + ! Deactivate coupling for active cycles? + call self % couplingInfo % moveToActive(self % activeTally) + call self % cycles(self % activeTally, self % activeAtch, self % N_active) + call self % couplingInfo % endCoupling() + ! Collect results from other processes call self % inactiveTally % collectDistributed() call self % activeTally % collectDistributed() @@ -334,6 +346,9 @@ subroutine cycles(self, tally, tallyAtch, N_cycles) call statusMsg("End time: " // trim(secToChar(end_T))) call statusMsg("Time to end: " // trim(secToChar(T_toEnd))) call tally % display() + + ! Perform coupling operations + call self % couplingInfo % couple(i) end do @@ -429,7 +444,7 @@ subroutine init(self, dict) type(visualiser) :: viz class(field), pointer :: field class(pieceConstantField), pointer :: pcField - real(defReal) :: maxDensityScale, maxTemperature + real(defReal) :: maxDensityScale, maxTemperature, maxCouplingTemp character(100), parameter :: Here ='init (eigenPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -449,6 +464,18 @@ subroutine init(self, dict) ! Parallel buffer size call dict % getOrDefault(self % bufferSize, 'buffer', 1000) + ! Check if the calculation is coupled + maxCouplingTemp = NO_TEMPERATURE + if (dict % isPresent('coupling')) then + tempDict => dict % getDictPtr('coupling') + call self % couplingInfo % init(tempDict) + call tempDict % getOrDefault(maxCouplingTemp, 'maxTemp', NO_TEMPERATURE) + if (.not. tempDict % isPresent('maxTemp')) then + call statusMsg('If coupling with temperature fields, it is advisable to include a maximum '& + //'possible temperature (maxTemp) in order to correctly precompute the majorant.') + end if + end if + ! Process type of data select case(energy) case('mg') @@ -537,6 +564,7 @@ subroutine init(self, dict) else maxTemperature = NO_TEMPERATURE end if + maxTemperature = max(maxTemperature, maxCouplingTemp) ! If present, build density field if (dict % isPresent('density')) then @@ -638,7 +666,11 @@ subroutine init(self, dict) ! Attach attachments to result tallies call self % inactiveTally % push(self % inactiveAtch) call self % activeTally % push(self % activeAtch) - + + ! Attach a tally admin for coupling + if (self % couplingInfo % doCoupling()) then + call self % couplingInfo % attachTally(self % inactiveTally) + end if call self % printSettings() diff --git a/SharedModules/timer_mod.f90 b/SharedModules/timer_mod.f90 index c87e57ba8..5735115cc 100644 --- a/SharedModules/timer_mod.f90 +++ b/SharedModules/timer_mod.f90 @@ -22,11 +22,14 @@ !! timerReset -> Reset a defined timer !! timerTime -> Get total elapsed time in Seconds !! +!! c_usleep -> pauses SCONE for several microseconds +!! !! module timer_mod use numPrecision use genericProcedures, only : fatalError + use iso_c_binding, only : c_int implicit none private @@ -38,6 +41,7 @@ module timer_mod public :: timerStop public :: timerReset public :: timerTime + public :: c_usleep !! !! This derived type allows to measure time in a program @@ -77,6 +81,16 @@ module timer_mod real(defReal),parameter :: GROWTH_RATIO = 1.6_defReal integer(shortInt),parameter :: MIN_SIZE = 5 + ! Interface to C usleep procedure + ! Takes microseconds as input + interface + subroutine c_usleep(useconds) bind(C, name="usleep") + import :: c_int + integer(c_int), value :: useconds + end subroutine + end interface + + contains !! @@ -338,4 +352,5 @@ function toSec_stopWatch(self) result(sec) end function toSec_stopWatch + end module timer_mod diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 69d92f068..9361bbbc8 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -18,11 +18,9 @@ module universalVariables integer(shortInt), parameter, public :: MAX_COL = 70 ! Maximum number of columns in console display ! Define variables which are important for tracking neutrons in the geometry - real(defReal), parameter, public :: INFINITY = 2.0_defReal**63, & - surface_tol = 1.0e-12_defReal, & ! Tol. on closeness to surface - SURF_TOL = 1.0E-12_defReal, & - INF = 2.0_defReal**63, & - NUDGE = 1.0e-8_defReal ! Distance to poke neutrons across boundaries for surface tracking + real(defReal), parameter, public :: INF = 2.0_defReal**63, & + SURF_TOL = 1.0E-12_defReal, & + NUDGE = 1.0e-8_defReal ! Distance to poke neutrons across boundaries for surface tracking ! Energy variables real(defReal), parameter, public :: MINIMUM_ENERGY = 1.0E-11_defReal, & diff --git a/Tallies/scoreMemory_class.f90 b/Tallies/scoreMemory_class.f90 index bf67940c5..c125ed1fc 100644 --- a/Tallies/scoreMemory_class.f90 +++ b/Tallies/scoreMemory_class.f90 @@ -66,6 +66,8 @@ module scoreMemory_class !! !! reduceBins(): Move the scores from parallelBins and different processes to bins. !! + !! reset(): Resets the memory: zeros all scores and reverts cycles/batches to 0. + !! !! Example use case: !! !! do batches=1,20 @@ -112,6 +114,7 @@ module scoreMemory_class procedure :: getBatchSize procedure :: reduceBins procedure :: collectDistributed + procedure :: reset ! Private procedures procedure, private :: score_defReal @@ -619,4 +622,17 @@ elemental function getScore(self, idx) result (score) end function getScore + !! + !! Resets score memory to having no scores, cycles, or batches + !! + subroutine reset(self) + class(scoreMemory), intent(inout) :: self + + self % bins = ZERO + self % parallelBins = ZERO + self % batchN = 0 + self % cycles = 0 + + end subroutine reset + end module scoreMemory_class diff --git a/Tallies/tallyAdmin_class.f90 b/Tallies/tallyAdmin_class.f90 index feba84a6d..839c5d111 100644 --- a/Tallies/tallyAdmin_class.f90 +++ b/Tallies/tallyAdmin_class.f90 @@ -67,17 +67,19 @@ module tallyAdmin_class !! push -> Push new attachment to the end of the list !! pop -> Remove attachment from the end of the list !! getEnd -> Returns pointer to the ned of the list (does not remove it like pop) - !! reportInColl -> Process pre-collision reports in all clerks - !! reportOutColl -> Process post-collision reports in all clerks - !! reportPath -> Process pathlength reports in all clerks - !! reportTrans -> Process transition reports in all clerks - !! reportSpawn -> Process particle generation reports in all clerks - !! reportHist -> Process history reports in all clerks - !! reportCycleStart -> Process start of cycle reports in all clerks - !! reportCycleEnd -> Process end of cycle reports in all clerks - !! getResult -> Return tallyResult object from a named Clerk + !! reportInColl -> Process pre-collision reports in all clerks + !! reportOutColl -> Process post-collision reports in all clerks + !! reportPath -> Process pathlength reports in all clerks + !! reportTrans -> Process transition reports in all clerks + !! reportSpawn -> Process particle generation reports in all clerks + !! reportHist -> Process history reports in all clerks + !! reportCycleStart -> Process start of cycle reports in all clerks + !! reportCycleEnd -> Process end of cycle reports in all clerks + !! getResult -> Return tallyResult object from a named Clerk + !! collectDistributed -> Sums tally results from all ranks into the master + !! resetMemory -> Resets score memory for all tallies !! display -> Call "display" on all clerks registered to display - !! isConverged -> Return .true. if all convergance targets have been reached + !! isConverged -> Return .true. if all convergence targets have been reached !! print -> Prints results to an output file object !! !! SAMPLE DICTIOANRY INPUT: @@ -151,6 +153,7 @@ module tallyAdmin_class ! Interaction procedures procedure :: getResult procedure :: collectDistributed + procedure :: resetMemory ! Display procedures procedure :: display @@ -410,7 +413,7 @@ end subroutine display !! None !! !! Result: - !! True is all convergance targets have been reached. False otherwise. + !! True if all convergence targets have been reached. False otherwise. !! !! Errors: !! None @@ -893,4 +896,14 @@ subroutine addToReports(self, reportCode, idx) end subroutine addToReports + !! + !! Resets memory to a fresh state + !! + subroutine resetMemory(self) + class(tallyAdmin),intent(inout) :: self + + call self % mem % reset() + + end subroutine resetMemory + end module tallyAdmin_class diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 8aefa39d7..345b50c63 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -91,12 +91,14 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle is in a region with overlapping cells case(OVERLAP_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index c1461b4a0..3434db458 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -131,11 +131,13 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle somehow ended in an overlap material case(OVERLAP_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default @@ -195,8 +197,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else @@ -251,12 +253,14 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle ended in an overlap material case(OVERLAP_MAT) - print*, 'Particle location: ', p % rGlobal() + print *, 'Particle location: ', p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index f73d27563..63055bfbf 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -70,8 +70,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else @@ -134,11 +134,13 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! Give error if the particle somehow ended in an undefined material case(UNDEF_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in undefined material") ! Give error if the particle is in a region with overlapping cells case(OVERLAP_MAT) print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() call fatalError(Here, "Particle is in overlapping cells") case default diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index 87285f94d..fb069e907 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -62,6 +62,9 @@ Example: :: uniformFissionSites { } varianceReduction { } source { } + coupling { } + temperature { } + density { } .. note:: Although a ``source`` definition is not required, it can be included to replace @@ -121,6 +124,8 @@ Example: :: *Optional entries* :: varianceReduction { } + temperature { } + density { } kineticPhysicsPackage ##################### @@ -190,6 +195,8 @@ Example: :: *Optional entries* :: varianceReduction { } + temperature { } + density { } alphaPhysicsPackage ################### @@ -249,6 +256,8 @@ Example: :: uniformFissionSites { } varianceReduction { } source { } + temperature { } + density { } .. note:: Although a ``source`` definition is not required, it can be included to replace @@ -987,6 +996,57 @@ Example: :: root { id 1000; type rootUniverse; border 10; fill u<1>; } +Coupling +-------- + +To performed coupled physics simulations (currently limited to k and alpha eigenvalue calculations), the keyword +``coupling`` must be present in the input file with an appropriately defined dictionary. The entries +in this dictionary are: +* tally: a tallyAdmin definition, used to produce tallies and normalisation required for coupling, e.g., power. +* receiveFile: tex file that SCONE can expect to receive from external codes. This will be read to signal either that + coupling is continuing (containing 'SIGUSR1') or that coupling is finished (containing 'SIGTERM'). + The latter could be based on some convergence criteria being met by the external solver. +* sendFile: text file that SCONE should write to send signals to other codes. This will signal whether coupling should + continue (containing 'SIGUSR1') or that coupling has finished (containing 'SIGTERM'). +* updateFreq: sets how often SCONE exchanges data with coupled codes. More frequently may mean faster convergence + of the coupled system, but results sent to other codes will be noisier due to poorer statistics. +* outputFile: name of the file containing tallies that SCONE will produce. +* outputFormat (*optional*, default = ``asciiMATLAB``): type of output file. + Choices are ``asciiMATLAB`` and ``asciiJSON``. +* fieldPaths: array of paths to files that SCONE will read which contain field definitions, produced by external solvers. + There must be as many paths as fields, and their ordering must match. +* fields: array of names of fields that SCONE will update during coupling. Can be ``temperature`` and/or ``density``. +* maxTemp: (*optional*, defaults to the highest material temperature initially present in the geometry): + the maximum temperature expected to be produced during coupling. This is important to correctly initialise the + majorant if using delta tracking. If this temperature is exceeded during coupling, a warning will be produced + regarding the validity of the results. +* maxIt: (*optional*, defaults to infinity) the maximum number of MC cycles that will take place before SCONE + terminates the coupling. SCONE will continue to run, but it will signal to the external solver that coupling + has finished and fields will no longer be updated. +* endActive: (*optional*, defaults to false) decides whether to end coupling during active cycles. This can save + computational time due to the cost of communication, but it may also produce higher variance/bias results. + +An example of a coupling dictionary is: :: + + coupling { + tally {norm power; normVal 100; + power {type collisionClerk; response (kappa) kappa {type macroResponse; MT -80;}} + powerLocal {type collisionClerk; response (kappa) kappa {type macroResponse; MT -80;} + map {type spaceMap; axis z; grid lin; min -200; min 200; N 40;} + } + } + receiveFile /path/toSCONE; + sendFile /path/fromSCONE; + updateFreq 5; + outputFile /path/tallies_from_SCONE; + outputFormat asciiMATLAB; + fieldPaths (/path1/tempField /path2/densityField); + fields (temperature density); + maxTemp 2000; + maxIt 100; + endActive 1; + } + Visualiser ----------