diff --git a/CollisionOperator/CollisionProcessors/CMakeLists.txt b/CollisionOperator/CollisionProcessors/CMakeLists.txt index ee923a0ed..8b8ab513a 100644 --- a/CollisionOperator/CollisionProcessors/CMakeLists.txt +++ b/CollisionOperator/CollisionProcessors/CMakeLists.txt @@ -2,5 +2,6 @@ add_sources( ./collisionProcessor_inter.f90 ./collisionProcessorFactory_func.f90 ./neutronCEstd_class.f90 - ./neutronCEimp_class.f90 - ./neutronMGstd_class.f90) + ./neutronCEimp_class.f90 + ./neutronMGstd_class.f90 + ./neutronMGimp_class.f90) diff --git a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 index 83692f38f..57ac0e188 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 @@ -11,6 +11,7 @@ module collisionProcessorFactory_func use neutronCEstd_class, only : neutronCEstd use neutronCEimp_class, only : neutronCEimp use neutronMGstd_class, only : neutronMGstd + use neutronMGimp_class, only : neutronMGimp implicit none private @@ -23,7 +24,8 @@ module collisionProcessorFactory_func ! For now it is necessary to adjust trailing blanks so all enteries have the same length character(nameLen),dimension(*),parameter :: AVALIBLE_collisionProcessors = [ 'neutronCEstd',& 'neutronCEimp',& - 'neutronMGstd'] + 'neutronMGstd',& + 'neutronMGimp'] contains @@ -54,6 +56,10 @@ subroutine new_collisionProcessor(new,dict) case('neutronMGstd') allocate(neutronMGstd :: new) + case('neutronMGimp') + allocate(neutronMGimp :: new) + call new % init(dict) + case default print *, AVALIBLE_collisionProcessors call fatalError(Here, 'Unrecognised type of collisionProcessor: ' // trim(type)) diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 index fabc11b3b..466951106 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 @@ -86,13 +86,15 @@ module collisionProcessor_inter !! Procedure interface for all customisable actions associated with !! processing of sollision event (scatter, fission etc.) !! - subroutine collisionAction(self, p, collDat, thisCycle, nextCycle) + subroutine collisionAction(self, p, tally, collDat, thisCycle, nextCycle) import :: collisionProcessor, & collisionData, & + tallyAdmin, & particle,& particleDungeon class(collisionProcessor), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -104,7 +106,7 @@ end subroutine collisionAction !! !! Generic flow of collision processing !! - subroutine collide(self, p, tally ,thisCycle, nextCycle) + subroutine collide(self, p, tally, thisCycle, nextCycle) class(collisionProcessor), intent(inout) :: self class(particle), intent(inout) :: p type(tallyAdmin), intent(inout) :: tally @@ -118,30 +120,32 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle) ! Report in-collision & save pre-collison state ! Note: the ordering must not be changed between feeding the particle to the tally - ! and updating the particle's preCollision state, otherwise this may cause certain + ! and updating the particle's preCollision state, otherwise this may cause certain ! tallies (e.g., collisionProbability) to return dubious results + call tally % reportInColl(p, .false.) + call p % savePreCollision() ! Choose collision nuclide and general type (Scatter, Capture or Fission) - call self % sampleCollision(p, collDat, thisCycle, nextCycle) + call self % sampleCollision(p, tally, collDat, thisCycle, nextCycle) ! Perform implicit treatment - call self % implicit(p, collDat, thisCycle, nextCycle) + call self % implicit(p, tally, collDat, thisCycle, nextCycle) ! Select physics to be processed based on MT number select case(collDat % MT) case(N_N_elastic, macroAllScatter) - call self % elastic(p, collDat, thisCycle, nextCycle) + call self % elastic(p, tally, collDat, thisCycle, nextCycle) case(N_N_inelastic, macroIEScatter) - call self % inelastic(p, collDat, thisCycle, nextCycle) + call self % inelastic(p, tally, collDat, thisCycle, nextCycle) case(N_DISAP, macroCapture) - call self % capture(p, collDat, thisCycle, nextCycle) + call self % capture(p, tally, collDat, thisCycle, nextCycle) case(N_FISSION, macroFission) - call self % fission(p, collDat, thisCycle, nextCycle) + call self % fission(p, tally, collDat, thisCycle, nextCycle) case(noInteraction) ! Do nothing @@ -152,7 +156,7 @@ subroutine collide(self, p, tally ,thisCycle, nextCycle) end select ! Apply post collision implicit treatments - call self % cutoffs(p, collDat, thisCycle, nextCycle) + call self % cutoffs(p, tally, collDat, thisCycle, nextCycle) ! Update particle collision counter p % collisionN = p % collisionN + 1 diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index 0d4f05acb..89d5a09d7 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 @@ -38,6 +38,10 @@ module neutronCEimp_class ! Scattering procedures use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, & asymptoticInelasticScatter, targetVelocity_DBRCXS + + ! Tally interfaces + use tallyAdmin_class, only : tallyAdmin + implicit none private @@ -58,6 +62,7 @@ module neutronCEimp_class !! impAbs -> is implicit capture performed? (off by default) !! impGen -> are fission sites generated implicitly? (on by default) !! UFS -> uniform fission sites variance reduction + !! maxSplit -> maximum number of splits allowed per particle (default = 1000) !! thresh_E -> Energy threshold for explicit treatment of target nuclide movement [-]. !! Target movment is sampled if neutron energy E < kT * thresh_E where !! kT is target material temperature in [MeV]. (default = 400.0) @@ -85,6 +90,7 @@ module neutronCEimp_class !! #impGen ;# !! #UFS ;# !! #weightWindows ;# + !! #maxSplit ;# !! } !! type, public, extends(collisionProcessor) :: neutronCEimp @@ -105,14 +111,15 @@ module neutronCEimp_class real(defReal) :: thresh_A real(defReal) :: DBRCeMin real(defReal) :: DBRCeMax + integer(shortInt) :: maxSplit ! Variance reduction options - logical(defBool) :: weightWindows - logical(defBool) :: splitting - logical(defBool) :: roulette - logical(defBool) :: implicitAbsorption ! Prevents particles dying through capture - logical(defBool) :: implicitSites ! Generates fission sites on every fissile collision - logical(defBool) :: uniFissSites + logical(defBool) :: weightWindows + logical(defBool) :: splitting + logical(defBool) :: roulette + logical(defBool) :: implicitAbsorption ! Prevents particles dying through capture + logical(defBool) :: implicitSites ! Generates fission sites on every fissile collision + logical(defBool) :: uniFissSites type(weightWindowsField), pointer :: weightWindowsMap @@ -164,6 +171,7 @@ subroutine init(self, dict) ! Obtain settings for variance reduction call dict % getOrDefault(self % weightWindows,'weightWindows', .false.) + call dict % getOrDefault(self % maxSplit,'maxSplit', 1000) call dict % getOrDefault(self % splitting,'split', .false.) call dict % getOrDefault(self % roulette,'roulette', .false.) call dict % getOrDefault(self % minWgt,'minWgt',0.25_defReal) @@ -190,8 +198,8 @@ subroutine init(self, dict) end if if (self % implicitAbsorption) then - if (.not.self % roulette) call fatalError(Here,& - 'Must use Russian roulette when using implicit absorption') + if (.not.self % roulette .and. .not. self % weightWindows) call fatalError(Here,& + 'Must use Russian roulette or weight windows when using implicit absorption') if (.not.self % implicitSites) call fatalError(Here,& 'Must generate fission sites implicitly when using implicit absorption') end if @@ -213,9 +221,10 @@ end subroutine init !! !! Samples collision without any implicit treatment !! - subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle) + subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -252,9 +261,10 @@ end subroutine sampleCollision !! !! Perform implicit treatment !! - subroutine implicit(self, p, collDat, thisCycle, nextCycle) + subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -302,7 +312,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle) ! Store new sites in the next cycle dungeon r = p % rGlobal() - do i=1,n + do i = 1,n call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG) dir = rotateVector(p % dirGlobal(), mu, phi) @@ -321,6 +331,9 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle) call nextCycle % detain(pTemp) if (self % uniFissSites) call self % ufsField % storeFS(pTemp) + ! Report birth of new particle + call tally % reportSpawn(N_FISSION, p, pTemp) + end do end if @@ -346,9 +359,10 @@ end subroutine implicit !! !! Process capture reaction !! - subroutine capture(self, p, collDat, thisCycle, nextCycle) + subroutine capture(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -360,9 +374,10 @@ end subroutine capture !! !! Process fission reaction !! - subroutine fission(self, p, collDat, thisCycle, nextCycle) + subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -426,6 +441,9 @@ subroutine fission(self, p, collDat, thisCycle, nextCycle) call nextCycle % detain(pTemp) if (self % uniFissSites) call self % ufsField % storeFS(pTemp) + ! Report birth of new particle + call tally % reportSpawn(N_FISSION, p, pTemp) + end do end if @@ -438,9 +456,10 @@ end subroutine fission !! !! All CE elastic scattering happens in the CM frame !! - subroutine elastic(self, p, collDat, thisCycle, nextCycle) + subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -476,9 +495,10 @@ end subroutine elastic !! !! Process inelastic scattering !! - subroutine inelastic(self, p, collDat, thisCycle, nextCycle) + subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -506,9 +526,10 @@ end subroutine inelastic !! !! Apply cutoffs !! - subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) + subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -516,7 +537,8 @@ subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) real(defReal) :: minWgt, maxWgt, avWgt if (p % isDead) then - ! Do nothing + + ! Do nothing ! elseif (p % E < self % minE) then p % isDead = .true. @@ -530,15 +552,18 @@ subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) ! If a particle is outside the WW map and all the weight limits ! are zero nothing happens. NOTE: this holds for positive weights only - if ((p % w > maxWgt) .and. (maxWgt /= ZERO)) then - call self % split(p, thisCycle, maxWgt) + + if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < self % maxSplit)) then + call self % split(p, tally, thisCycle, maxWgt) + elseif (p % w < minWgt) then call self % russianRoulette(p, avWgt) + end if ! Splitting with fixed threshold elseif ((self % splitting) .and. (p % w > self % maxWgt)) then - call self % split(p, thisCycle, self % maxWgt) + call self % split(p, tally, thisCycle, self % maxWgt) ! Roulette with fixed threshold and survival weight elseif ((self % roulette) .and. (p % w < self % minWgt)) then @@ -567,24 +592,41 @@ end subroutine russianRoulette !! !! Split particle which has too large a weight !! - subroutine split(self, p, thisCycle, maxWgt) + subroutine split(self, p, tally, thisCycle, maxWgt) class(neutronCEimp), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally class(particleDungeon), intent(inout) :: thisCycle real(defReal), intent(in) :: maxWgt + type(particleState) :: pTemp integer(shortInt) :: mult, i ! This value must be at least 2 mult = ceiling(p % w/maxWgt) - ! Decrease weight - p % w = p % w/mult + ! Limit maximum split + if (mult + p % splitCount > self % maxSplit) then + mult = self % maxSplit - p % splitCount + 1 + end if + + ! Copy particle to a particle state + ! Note that particleState doesn't have property splitCount, so it is reset + ! to 0 for the new particle + pTemp = p + pTemp % wgt = p % w/mult ! Add split particle's to the dungeon - do i = 1,mult-1 - call thisCycle % detain(p) + do i = 1, mult-1 + call thisCycle % detain(pTemp) + call tally % reportSpawn(N_N_SPLIT, p, pTemp) end do + ! Update particle split count + p % splitCount = p % splitCount + mult + + ! Decrease original particle weight + p % w = p % w/mult + end subroutine split !! diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 index 7070a456f..e2ddf2f44 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 @@ -32,6 +32,10 @@ module neutronCEstd_class ! Scattering procedures use scatteringKernels_func, only : asymptoticScatter, targetVelocity_constXS, & asymptoticInelasticScatter, targetVelocity_DBRCXS + + ! Tally interfaces + use tallyAdmin_class, only : tallyAdmin + implicit none private @@ -135,9 +139,10 @@ end subroutine init !! !! Samples collision without any implicit treatment !! - subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle) + subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -174,9 +179,10 @@ end subroutine sampleCollision !! !! Perform implicit treatment !! - subroutine implicit(self, p, collDat, thisCycle, nextCycle) + subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -216,7 +222,7 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle) wgt = sign(w0, wgt) r = p % rGlobal() - do i=1,n + do i = 1,n call fission % sampleOut(mu, phi, E_out, p % E, p % pRNG) dir = rotateVector(p % dirGlobal(), mu, phi) @@ -233,6 +239,10 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle) pTemp % collisionN = 0 call nextCycle % detain(pTemp) + + ! Report birth of new particle + call tally % reportSpawn(N_FISSION, p, pTemp) + end do end if @@ -241,9 +251,10 @@ end subroutine implicit !! !! Process capture reaction !! - subroutine capture(self, p, collDat, thisCycle, nextCycle) + subroutine capture(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -255,9 +266,10 @@ end subroutine capture !! !! Process fission reaction !! - subroutine fission(self, p, collDat, thisCycle, nextCycle) + subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -271,9 +283,10 @@ end subroutine fission !! !! All CE elastic scattering happens in the CM frame !! - subroutine elastic(self, p, collDat, thisCycle, nextCycle) + subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -309,9 +322,10 @@ end subroutine elastic !! !! Process inelastic scattering !! - subroutine inelastic(self, p, collDat, thisCycle, nextCycle) + subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -339,9 +353,10 @@ end subroutine inelastic !! !! Apply cutoffs !! - subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) + subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle) class(neutronCEstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle diff --git a/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 new file mode 100644 index 000000000..bc2155c22 --- /dev/null +++ b/CollisionOperator/CollisionProcessors/neutronMGimp_class.f90 @@ -0,0 +1,393 @@ +module neutronMGimp_class + + use numPrecision + use endfConstants + use universalVariables, only : nameWW + use genericProcedures, only : fatalError, rotateVector, numToChar + use dictionary_class, only : dictionary + use RNG_class, only : RNG + + ! Particle types + use particle_class, only : particle, particleState, printType, P_NEUTRON + use particleDungeon_class, only : particleDungeon + + ! Abstract interface + use collisionProcessor_inter, only : collisionProcessor, collisionData ,init_super => init + + ! Nuclear Data Interface + use nuclearDataReg_mod, only : ndReg_getNeutronMG => getNeutronMG + use nuclearDatabase_inter, only : nuclearDatabase + use mgNeutronDatabase_inter, only : mgNeutronDatabase + use mgNeutronMaterial_inter, only : mgNeutronMaterial, mgNeutronMaterial_CptrCast + use reactionHandle_inter, only : reactionHandle + use multiScatterMG_class, only : multiScatterMG, multiScatterMG_CptrCast + use fissionMG_class, only : fissionMG, fissionMG_TptrCast + + ! Cross section packages + use neutronXsPackages_class, only : neutronMacroXSs + + ! Geometry and fields + use geometryReg_mod, only : gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + use weightWindowsField_class, only : weightWindowsField, weightWindowsField_TptrCast + + ! Tally interfaces + use tallyAdmin_class, only : tallyAdmin + + implicit none + private + + !! + !! Scalar collision processor for MG neutrons + !! -> Preforms implicit fission site generation + !! -> Preforms analog capture + !! -> Treats fission as capture (only implicit generation of 2nd-ary neutrons) + !! -> Does not create secondary non-neutron projectiles + !! -> Supports the use of weight windows + !! + !! Settings: + !! weightWindows -> uses a weight windows field (off by default) + !! maxSplit -> maximum number of splits allowed per particle (default = 1000) + !! + !! Sample dictionary input: + !! collProcName { + !! type neutronMGimp; + !! #weightWindows ;# + !! #maxSplit ;# + !! } + !! + type, public, extends(collisionProcessor) :: neutronMGimp + private + class(mgNeutronDatabase), pointer, public :: xsData => null() + class(mgNeutronMaterial), pointer, public :: mat => null() + + !! Settings - private + integer(shortInt) :: maxSplit + + ! Variance reduction options + logical(defBool) :: weightWindows + type(weightWindowsField), pointer :: weightWindowsMap + + contains + ! Initialisation procedure + procedure :: init + + ! Implementation of customisable procedures + procedure :: sampleCollision + procedure :: implicit + procedure :: elastic + procedure :: inelastic + procedure :: capture + procedure :: fission + procedure :: cutoffs + + ! Variance reduction procedures + procedure, private :: split + procedure, private :: russianRoulette + + end type neutronMGimp + +contains + + !! + !! Initialise from dictionary + !! + subroutine init(self, dict) + class(neutronMGimp), intent(inout) :: self + class(dictionary), intent(in) :: dict + integer(shortInt) :: idx + character(100), parameter :: Here = 'init (neutronMGimp_class.f90)' + + ! Call superclass + call init_super(self, dict) + + ! Obtain settings for variance reduction + call dict % getOrDefault(self % maxSplit,'maxSplit', 1000) + call dict % getOrDefault(self % weightWindows,'weightWindows', .false.) + + ! Sets up the weight windows field + if (self % weightWindows) then + idx = gr_fieldIdx(nameWW) + self % weightWindowsMap => weightWindowsField_TptrCast(gr_fieldPtr(idx)) + end if + + end subroutine init + + !! + !! Samples collision without any implicit treatment + !! + subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + type(neutronMacroXSs) :: macroXSs + real(defReal) :: r + character(100),parameter :: Here =' sampleCollision (neutronMGimp_class.f90)' + + ! Verify that particle is MG neutron + if( .not. p % isMG .or. p % type /= P_NEUTRON) then + call fatalError(Here, 'Supports only MG Neutron. Was given CE '//printType(p % type)) + end if + + ! Verify and load nuclear data pointer + self % xsData => ndReg_getNeutronMG() + if(.not.associated(self % xsData)) call fatalError(Here, "Failed to get active database for MG Neutron") + + ! Get and verify material pointer + self % mat => mgNeutronMaterial_CptrCast( self % xsData % getMaterial( p % matIdx())) + if(.not.associated(self % mat)) call fatalError(Here, "Failed to get MG Neutron Material") + + ! Select Main reaction channel + call self % mat % getMacroXSs(macroXSs, p % G, p % pRNG) + r = p % pRNG % get() + + collDat % MT = macroXSs % invert(r) + + end subroutine sampleCollision + + !! + !! Preform implicit treatment + !! + subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + type(neutronMacroXSs) :: macroXSs + type(fissionMG),pointer :: fission + type(particleState) :: pTemp + real(defReal),dimension(3) :: r, dir + integer(shortInt) :: G_out, n, i + real(defReal) :: wgt, w0, rand1, mu, phi + real(defReal) :: sig_tot, k_eff, sig_nufiss + character(100),parameter :: Here = 'implicit (neutronMGimp_class.f90)' + + if ( self % mat % isFissile()) then + ! Obtain required data + wgt = p % w ! Current weight + w0 = p % preHistory % wgt ! Starting weight + k_eff = p % k_eff ! k_eff for normalisation + rand1 = p % pRNG % get() ! Random number to sample sites + + call self % mat % getMacroXSs(macroXSs, p % G, p % pRNG) + + sig_tot = macroXSs % total + sig_nuFiss = macroXSs % nuFission + + ! Sample number of fission sites generated + !n = int(wgt * sig_nuFiss/(sig_tot*k_eff) + r1, shortInt) + n = int(abs( (wgt * sig_nuFiss) / (w0 * sig_tot * k_eff)) + rand1, shortInt) + + ! Shortcut if no particles were samples + if (n < 1) return + + ! Get Fission reaction object + fission => fissionMG_TptrCast( self % xsData % getReaction(macroFission, collDat % matIdx)) + if (.not.associated(fission)) call fatalError(Here, 'Failed to getrive fissionMG reaction object') + + ! Store new sites in the next cycle dungeon + wgt = sign(w0, wgt) + r = p % rGlobal() + + do i=1,n + call fission % sampleOut(mu, phi, G_out, p % G, p % pRNG) + dir = rotateVector(p % dirGlobal(), mu, phi) + + ! Copy extra detail from parent particle (i.e. time, flags ect.) + pTemp = p + + ! Overwrite position, direction, energy group and weight + pTemp % r = r + pTemp % dir = dir + pTemp % G = G_out + pTemp % wgt = wgt + + call nextCycle % detain(pTemp) + + ! Report birth of new particle + call tally % reportSpawn(N_FISSION, p, pTemp) + + end do + end if + + end subroutine implicit + + !! + !! Elastic Scattering + !! + subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + ! Do nothing. Should not be called + + end subroutine elastic + + !! + !! Preform scattering + !! + subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + class(multiScatterMG),pointer :: scatter + integer(shortInt) :: G_out ! Post-collision energy group + real(defReal) :: phi ! Azimuthal scatter angle + real(defReal) :: w_mul ! Weight multiplier + character(100),parameter :: Here = "inelastic (neutronMGimp_class.f90)" + + ! Assign MT number + collDat % MT = macroIEscatter + + ! Get Scatter object + scatter => multiScatterMG_CptrCast( self % xsData % getReaction(macroIEscatter, collDat % matIdx)) + if(.not.associated(scatter)) call fatalError(Here, "Failed to get scattering reaction object for MG neutron") + + ! Sample Mu and G_out + call scatter % sampleOut(collDat % muL, phi, G_out, p % G, p % pRNG) + + ! Read scattering multiplicity + w_mul = scatter % production(p % G, G_out) + + ! Update neutron state + p % G = G_out + p % w = p % w * w_mul + call p % rotate(collDat % muL, phi) + + end subroutine inelastic + + !! + !! Preform capture + !! + subroutine capture(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + p % isDead = .true. + + end subroutine capture + + !! + !! Preform fission + !! + subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + p % isDead = .true. + + end subroutine fission + + !! + !! Applay cutoffs or post-collision implicit treatment + !! + subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + real(defReal), dimension(3) :: val + real(defReal) :: minWgt, maxWgt, avWgt + + if (p % isDead) then + ! Do nothing ! + + ! Weight Windows treatment + elseif (self % weightWindows) then + val = self % weightWindowsMap % at(p) + minWgt = val(1) + maxWgt = val(2) + avWgt = val(3) + + ! If a particle is outside the WW map and all the weight limits + ! are zero nothing happens. NOTE: this holds for positive weights only + if ((p % w > maxWgt) .and. (maxWgt /= ZERO) .and. (p % splitCount < self % maxSplit)) then + call self % split(p, tally, thisCycle, maxWgt) + + elseif (p % w < minWgt) then + call self % russianRoulette(p, avWgt) + + end if + + end if + + end subroutine cutoffs + + !! + !! Perform Russian roulette on a particle + !! + subroutine russianRoulette(self, p, avWgt) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal), intent(in) :: avWgt + + if (p % pRNG % get() < (ONE - p % w/avWgt)) then + p % isDead = .true. + else + p % w = avWgt + end if + + end subroutine russianRoulette + + !! + !! Split particle which has too large a weight + !! + subroutine split(self, p, tally, thisCycle, maxWgt) + class(neutronMGimp), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + class(particleDungeon), intent(inout) :: thisCycle + real(defReal), intent(in) :: maxWgt + type(particleState) :: pTemp + integer(shortInt) :: mult, i + + ! This value must be at least 2 + mult = ceiling(p % w/maxWgt) + + ! Limit maximum split + if (mult + p % splitCount > self % maxSplit) then + mult = self % maxSplit - p % splitCount + 1 + end if + + ! Copy particle to a particle state + ! Note that particleState doesn't have property splitCount, so it is reset + ! to 0 for the new particle + pTemp = p + pTemp % wgt = p % w/mult + + ! Add split particle's to the dungeon + do i = 1, mult-1 + call thisCycle % detain(pTemp) + call tally % reportSpawn(N_N_SPLIT, p, pTemp) + end do + + ! Update particle split count + p % splitCount = p % splitCount + mult + + ! Decrease original particle weight + p % w = p % w/mult + + end subroutine split + +end module neutronMGimp_class diff --git a/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 index 0665c7891..b816c6fe6 100644 --- a/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronMGstd_class.f90 @@ -25,13 +25,8 @@ module neutronMGstd_class ! Cross section packages use neutronXsPackages_class, only : neutronMacroXSs - - ! Nuclear Data - !use nuclearData_inter, only : nuclearData - !use perMaterialNuclearDataMG_inter, only : perMaterialNuclearDataMG - - ! Cross-section packages to interface with nuclear data - !use xsMacroSet_class, only : xsMacroSet, xsMacroSet_ptr + ! Tally interfaces + use tallyAdmin_class, only : tallyAdmin implicit none private @@ -87,9 +82,10 @@ end subroutine init !! !! Samples collision without any implicit treatment !! - subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle) + subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -121,9 +117,10 @@ end subroutine sampleCollision !! !! Preform implicit treatment !! - subroutine implicit(self, p, collDat, thisCycle, nextCycle) + subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -178,6 +175,10 @@ subroutine implicit(self, p, collDat, thisCycle, nextCycle) pTemp % collisionN = 0 call nextCycle % detain(pTemp) + + ! Report birth of new particle + call tally % reportSpawn(N_FISSION, p, pTemp) + end do end if @@ -186,9 +187,10 @@ end subroutine implicit !! !! Elastic Scattering !! - subroutine elastic(self, p , collDat, thisCycle, nextCycle) + subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -200,9 +202,10 @@ end subroutine elastic !! !! Preform scattering !! - subroutine inelastic(self, p, collDat, thisCycle, nextCycle) + subroutine inelastic(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -235,9 +238,10 @@ end subroutine inelastic !! !! Preform capture !! - subroutine capture(self, p, collDat, thisCycle, nextCycle) + subroutine capture(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -249,9 +253,10 @@ end subroutine capture !! !! Preform fission !! - subroutine fission(self, p, collDat, thisCycle, nextCycle) + subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle @@ -263,9 +268,10 @@ end subroutine fission !! !! Applay cutoffs or post-collision implicit treatment !! - subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) + subroutine cutoffs(self, p, tally, collDat, thisCycle, nextCycle) class(neutronMGstd), intent(inout) :: self class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally type(collisionData), intent(inout) :: collDat class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle diff --git a/InputFiles/c5g7_mox b/InputFiles/c5g7_mox index b8e86a4af..dbb696fc8 100644 --- a/InputFiles/c5g7_mox +++ b/InputFiles/c5g7_mox @@ -106,6 +106,7 @@ viz { centre (0.0 0.0 0.0); //width (25.0 25.0); axis z; + offset -17; res (500 500); } } @@ -119,6 +120,7 @@ nuclearData { water { temp 75675; + rgb (0 0 139); // Colour of water is dark blue moder { 1001.03 h-h2o.42; } composition { 1001.03 6.700E-002; diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 21ac73d11..78280d7d0 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 diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 3236083d2..729e209da 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -1,13 +1,14 @@ !! -!! Material Menu is a module (Singleton) that contains global definitions of diffrent materials +!! Material Menu is a module (Singleton) that contains global definitions of different materials !! !! It exists to make it easier for all databases to refer to the same materials by the same -!! name and index. This is necessary to avoid confusion resulting from diffrent materials with the -!! same name or index in diffrent databases. +!! name and index. This is necessary to avoid confusion resulting from different materials with the +!! same name or index in different databases. !! !! Public Members: !! materialDefs -> array of material definitions of type materialItem !! nameMap -> Map that maps material name to matIdx +!! colourMap -> Map that maps matIdx to 24bit colour (to use for visualisation) !! !! Interface: !! init -> Load material definitions from a dictionary @@ -21,8 +22,10 @@ module materialMenu_mod use numPrecision - use universalVariables, only : NOT_FOUND, VOID_MAT, OUTSIDE_MAT + use universalVariables, only : NOT_FOUND, VOID_MAT, OUTSIDE_MAT, UNDEF_MAT use genericProcedures, only : fatalError, charToInt, numToChar + use colours_func, only : rgb24bit + use intMap_class, only : intMap use charMap_class, only : charMap use dictionary_class, only : dictionary @@ -33,9 +36,9 @@ module materialMenu_mod !! Information about a single nuclide !! !! Based somewhat on MCNP conventions. - !! Atomic and Mass number identify cleary a nuclide species + !! Atomic and Mass number identify clearly a nuclide species !! Evaluation number T allows to refer to multiple states/evaluations of the same nuclide species - !! E.G. at a Diffrent temperature as in MCNP Library. + !! E.G. at a Different temperature as in MCNP Library. !! !! Public members: !! Z -> Atomic number @@ -83,6 +86,7 @@ module materialMenu_mod !! 5010.03 2.0E-005; !! } !! xsFile /home/uberMoffTarkin/XS/mat1.xs; + !! #rgb (255 0 0); # // RGB colour to be used in visualisation !! } !! !! NOTE: the moder dictionary is optional, necessary only if S(a,b) thermal scattering @@ -102,9 +106,16 @@ module materialMenu_mod procedure :: display => display_materialItem end type materialItem -!! MODULE COMPONENTS + !! Parameters + integer(shortInt), parameter :: COL_OUTSIDE = int(z'ffffff', shortInt) + integer(shortInt), parameter :: COL_VOID = int(z'000000', shortInt) + integer(shortInt), parameter :: COL_UNDEF = int(z'00ff00', shortInt) + + + !! MODULE COMPONENTS type(materialItem),dimension(:),allocatable,target,public :: materialDefs - type(charMap),target,public :: nameMap + type(charMap), target, public :: nameMap + type(intMap), public :: colourMap public :: init public :: kill @@ -131,7 +142,7 @@ subroutine init(dict) integer(shortInt) :: i character(nameLen) :: temp - ! Clean whatever may be alrady present + ! Clean whatever may be already present call kill() ! Load all material names @@ -142,17 +153,20 @@ subroutine init(dict) ! Load definitions do i=1,size(matNames) - call materialDefs(i) % init(matNames(i), dict % getDictPtr(matNames(i))) - materialDefs(i) % matIdx = i + call materialDefs(i) % init(matNames(i), i, dict % getDictPtr(matNames(i))) call nameMap % add(matNames(i), i) end do - ! Add special Material keywords to thedictionary + ! Add special Material keywords to the dictionary temp = 'void' call nameMap % add(temp, VOID_MAT) temp = 'outside' call nameMap % add(temp, OUTSIDE_MAT) + !! Load colours for the special materials + call colourMap % add(VOID_MAT, COL_VOID) + call colourMap % add(OUTSIDE_MAT, COL_OUTSIDE) + call colourMap % add(UNDEF_MAT, COL_UNDEF) end subroutine init @@ -204,7 +218,7 @@ end subroutine display !! Result: !! nameLen long character with material name !! - !! Erorrs: + !! Error: !! If idx is -ve or larger then number of defined materials !! Empty string '' is returned as its name !! @@ -250,25 +264,30 @@ end function matIdx !! !! Args: !! name [in] -> character with material name + !! idx [in] -> material index !! dict [in] -> dictionary with material definition !! !! Errors: !! FatalError if dictionary does not contain valid material definition. !! - subroutine init_materialItem(self, name, dict) - class(materialItem), intent(inout) :: self - character(nameLen),intent(in) :: name - class(dictionary), intent(in) :: dict - character(nameLen),dimension(:),allocatable :: keys, moderKeys - integer(shortInt) :: i - class(dictionary),pointer :: compDict, moderDict - logical(defBool) :: hasSab + subroutine init_materialItem(self, name, idx, dict) + class(materialItem), intent(inout) :: self + character(nameLen), intent(in) :: name + integer(shortInt), intent(in) :: idx + class(dictionary), intent(in) :: dict + character(nameLen), dimension(:), allocatable :: keys, moderKeys + integer(shortInt), dimension(:), allocatable :: temp + integer(shortInt) :: i + class(dictionary),pointer :: compDict, moderDict + logical(defBool) :: hasSab + character(100), parameter :: Here = 'init_materialItem (materialMenu_mod.f90)' ! Return to initial state call self % kill() ! Load easy components c self % name = name + self % matIdx = idx call dict % get(self % T,'temp') ! Get composition dictionary and load composition @@ -300,6 +319,17 @@ subroutine init_materialItem(self, name, dict) call self % nuclides(i) % init(keys(i)) end do + ! Add colour info if present + if(dict % isPresent('rgb')) then + call dict % get(temp, 'rgb') + + if (size(temp) /= 3) then + call fatalError(Here, "'rgb' keyword must have 3 values") + end if + + call colourMap % add(idx, rgb24bit(temp(1), temp(2), temp(3))) + end if + ! Save dictionary self % extraInfo = dict @@ -392,7 +422,7 @@ function isNucDefinition(key) result(isIt) za = verify(key(1:L),'0123456789') tt = verify(key(1:L),'0123456789', back = .true.) - ! Verify that the location of the dot is consistant + ! Verify that the location of the dot is consistent isIt = dot == za .and. dot == tt end function isNucDefinition @@ -437,7 +467,7 @@ end subroutine init_nuclideInfo !! None !! !! Result: - !! Character in format ZZAAA.TT that dscribes nuclide definition + !! Character in format ZZAAA.TT that describes nuclide definition !! !! Errors: !! None @@ -461,17 +491,17 @@ end function toChar_nuclideInfo !! Get pointer to a material definition under matIdx !! !! Args: - !! matIdx [in] -> Index of the material + !! idx [in] -> Index of the material !! !! Result: !! Pointer to a materialItem with the definition !! !! Errors: - !! FatalError if matIdx does not correspond to any defined material + !! FatalError if idx does not correspond to any defined material !! FatalError if material definitions were not loaded !! - function getMatPtr(matIdx) result(ptr) - integer(shortInt), intent(in) :: matIdx + function getMatPtr(idx) result(ptr) + integer(shortInt), intent(in) :: idx type(materialItem), pointer :: ptr character(100), parameter :: Here = 'getMatPtr (materialMenu_mod.f90)' @@ -481,13 +511,13 @@ function getMatPtr(matIdx) result(ptr) end if ! Verify matIdx - if( matIdx <= 0 .or. matIdx > nMat()) then - call fatalError(Here,"matIdx: "//numToChar(matIdx)// & + if( idx <= 0 .or. idx > nMat()) then + call fatalError(Here,"matIdx: "//numToChar(idx)// & " does not correspond to any defined material") end if ! Attach pointer - ptr => materialDefs(matIdx) + ptr => materialDefs(idx) end function getMatPtr diff --git a/ParticleObjects/particleDungeon_class.f90 b/ParticleObjects/particleDungeon_class.f90 index 7d47337ec..4d77538dd 100644 --- a/ParticleObjects/particleDungeon_class.f90 +++ b/ParticleObjects/particleDungeon_class.f90 @@ -27,8 +27,8 @@ module particleDungeon_class !! Dungeon can work like stacks or arrays. Stack-like behaviour is not really thread safe !! so it can be utilised when collecting and processing secondary particles in history !! that should be processed during the course of one cycle. Alternatively, one can use the - !! critical variations of the stack-like procedures. - !! Array-like behaviour allows to easily distribute particles among threads. As long as indices + !! critical variations of the stack-like procedures. + !! Array-like behaviour allows to easily distribute particles among threads. As long as indices !! assigned to different threads do not overlap, reading is thread-safe (I hope-MAK). !! !! @@ -170,18 +170,18 @@ subroutine detainCritical_particle(self,p) ! Increase population and weight self % pop = self % pop + 1 pop = self % pop - + ! Check for population overflow if (pop > size(self % prisoners)) then call fatalError(Here,'Run out of space for particles.& & Max size:'//numToChar(size(self % prisoners)) //& ' Current population: ' // numToChar(self % pop)) end if - + ! Load new particle self % prisoners(pop) = p !$omp end critical (dungeon) - + end subroutine detainCritical_particle !! @@ -224,16 +224,16 @@ subroutine detainCritical_particleState(self,p_state) !$omp critical (dungeon) self % pop = self % pop + 1 pop = self % pop - + ! Check for population overflow if (pop > size(self % prisoners)) then call fatalError(Here,'Run out of space for particles.& & Max size:'//numToChar(size(self % prisoners)) //& ' Current population: ' // numToChar(self % pop)) end if - + ! Load new particle - self % prisoners(pop) = p_state + self % prisoners(pop) = p_state !$omp end critical (dungeon) end subroutine detainCritical_particleState @@ -272,11 +272,11 @@ subroutine releaseCritical(self, p) ! Decrease population pop = self % pop self % pop = self % pop - 1 - + ! Load data into the particle p = self % prisoners(pop) !$omp end critical (dungeon) - + p % isDead = .false. end subroutine releaseCritical @@ -396,18 +396,19 @@ end subroutine normWeight !! Randomly duplicate or remove particles to match the number. !! Does not take weight of a particle into account! !! - subroutine normSize(self,N,rand) + subroutine normSize(self, N, rand) class(particleDungeon), intent(inout) :: self integer(shortInt), intent(in) :: N class(RNG), intent(inout) :: rand - integer(shortInt) :: excessP + integer(shortInt) :: excessP, n_copies, n_duplicates integer(shortInt) :: i, idx + integer(shortInt), dimension(:), allocatable :: duplicates character(100), parameter :: Here =' normSize (particleDungeon_class.f90)' ! Protect against invalid N if( N > size(self % prisoners)) then call fatalError(Here,'Requested size: '//numToChar(N) //& - 'is greather then max size: '//numToChar(size(self % prisoners))) + 'is greater then max size: '//numToChar(size(self % prisoners))) else if ( N <= 0 ) then call fatalError(Here,'Requested size: '//numToChar(N) //' is not +ve') end if @@ -416,9 +417,9 @@ subroutine normSize(self,N,rand) excessP = self % pop - N if (excessP > 0 ) then ! Reduce population with reservoir sampling - do i=N,self % pop + do i = N + 1, self % pop ! Select new index. Copy data if it is in the safe zone (<= N). - idx = int(i * rand % get())+1 + idx = int(i * rand % get()) + 1 if (idx <= N) then self % prisoners(idx) = self % prisoners(i) end if @@ -426,9 +427,28 @@ subroutine normSize(self,N,rand) self % pop = N else if (excessP < 0) then ! Clone randomly selected particles - do i = self % pop, N - idx = int(self % pop * rand % get()) + 1 - self % prisoners(i) = self % prisoners(idx) + ! For a massive undersampling duplicate (or n-plicate) particles + excessP = -excessP + n_copies = excessP / self % pop + n_duplicates = modulo(excessP, self % pop) + + ! Copy all particle maximum possible number of times + do i = 1, n_copies + self % prisoners(N * i + 1 : N * (i + 1)) = self % prisoners(1:N) + end do + + ! Choose the remainder particles to duplicate without replacement + duplicates = [(i, i = 1, n_duplicates)] + do i = n_duplicates + 1, self % pop + idx = int(i * rand % get()) + 1 + if (idx <= n_duplicates) then + duplicates(idx) = i + end if + end do + + ! Copy the duplicated particles at the end + do i = 1, n_duplicates + self % prisoners(self % pop + i) = self % prisoners(duplicates(i)) end do self % pop = N diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index 0fde0fd91..69021c389 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -116,6 +116,7 @@ module particle_class class(RNG), pointer :: pRNG => null() ! Pointer to RNG associated with the particle real(defReal) :: k_eff ! Value of default keff for implicit source generation integer(shortInt) :: geomIdx ! Index of the geometry used by the particle + integer(shortInt) :: splitCount = 0 ! Counter of number of splits ! Archived snapshots of previous states type(particleState) :: preHistory @@ -273,6 +274,7 @@ pure subroutine particle_fromParticleState(LHS,RHS) LHS % type = RHS % type LHS % time = RHS % time LHS % collisionN = RHS % collisionN + LHS % splitCount = 0 ! Reinitialise counter for number of splits end subroutine particle_fromParticleState diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index aa89111c1..d67603e8a 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -323,7 +323,7 @@ subroutine collectResults(self) type(outputFile) :: out character(nameLen) :: name - call out % init(self % outputFormat) + call out % init(self % outputFormat, filename=self % outputFile) name = 'seed' call out % printValue(self % pRNG % getSeed(),name) @@ -359,8 +359,6 @@ subroutine collectResults(self) call self % activeTally % print(out) call out % endBlock() - call out % writeToFile(self % outputFile) - end subroutine collectResults diff --git a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 index 35e2ff2d3..b1dbb4e86 100644 --- a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 +++ b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 @@ -269,7 +269,7 @@ subroutine collectResults(self) type(outputFile) :: out character(nameLen) :: name - call out % init(self % outputFormat) + call out % init(self % outputFormat, filename=self % outputFile) name = 'seed' call out % printValue(self % pRNG % getSeed(),name) @@ -290,8 +290,6 @@ subroutine collectResults(self) ! Print tally call self % tally % print(out) - call out % writeToFile(self % outputFile) - end subroutine collectResults diff --git a/SharedModules/CMakeLists.txt b/SharedModules/CMakeLists.txt index dcd1633ec..2a8baa4ce 100644 --- a/SharedModules/CMakeLists.txt +++ b/SharedModules/CMakeLists.txt @@ -12,7 +12,8 @@ add_sources( ./genericProcedures.f90 ./timer_mod.f90 ./charLib_func.f90 ./openmp_func.f90 - ./errors_mod.f90) + ./errors_mod.f90 + ./colours_func.f90) add_unit_tests( ./Tests/grid_test.f90 ./Tests/energyGrid_test.f90 @@ -21,4 +22,5 @@ add_unit_tests( ./Tests/grid_test.f90 ./Tests/hashFunctions_test.f90 ./Tests/timer_test.f90 ./Tests/conversions_test.f90 - ./Tests/charLib_test.f90) + ./Tests/charLib_test.f90 + ./Tests/colours_test.f90) diff --git a/SharedModules/Tests/colours_test.f90 b/SharedModules/Tests/colours_test.f90 new file mode 100644 index 000000000..3aba28b02 --- /dev/null +++ b/SharedModules/Tests/colours_test.f90 @@ -0,0 +1,26 @@ +module colours_test + use numPrecision + use colours_func, only : rgb24bit + use pFUnit_mod + + implicit none + +contains + + +@Test + subroutine testColourConversions() + + ! Test by comparison with some hex values + @assertEqual(int(z"123456"), rgb24bit(18, 52, 86)) + + @assertEqual(int(z"000000"), rgb24bit(0, 0, 0)) + @assertEqual(int(z"ffffff"), rgb24bit(255, 255, 255)) + + @assertEqual(int(z"ff0000"), rgb24bit(255, 0, 0)) + @assertEqual(int(z"00ff00"), rgb24bit(0, 255, 0)) + @assertEqual(int(z"0000ff"), rgb24bit(0, 0, 255)) + + end subroutine testColourConversions + +end module colours_test diff --git a/SharedModules/colours_func.f90 b/SharedModules/colours_func.f90 new file mode 100644 index 000000000..5e7fec61a --- /dev/null +++ b/SharedModules/colours_func.f90 @@ -0,0 +1,48 @@ +!! +!! Contains function to deal with colours +!! +module colours_func + use numPrecision + use genericProcedures, only: numToChar + use errors_mod, only: fatalError + implicit none + + public :: rgb24bit + +contains + + + !! + !! Return integer with 24bit colour value from r,g,b components + !! + !! Args: + !! r [in] -> red component in [0-255] range + !! g [in] -> green component in [0-255] range + !! b [in] -> blue component in [0-255] range + !! + function rgb24bit(r, g, b) result(col) + integer(shortInt), intent(in) :: r + integer(shortInt), intent(in) :: g + integer(shortInt), intent(in) :: b + integer(shortInt) :: col + character(100), parameter :: Here = "rgb24bit (colours_func.f90)" + + ! Check that the values are in the correct range + if (r < 0 .or. r > 255) then + call fatalError(Here, "Red value is out of [0-255] range: " // numToChar(r)) + end if + if (g < 0 .or. g > 255) then + call fatalError(Here, "Green value is out of [0-255] range: " // numToChar(g)) + end if + if (b < 0 .or. b > 255) then + call fatalError(Here, "Blue value is out of [0-255] range: " // numToChar(b)) + end if + + ! Compute the 24 bit colour + ! Note inversion, blue is the least significant bits + col = b + 256 * g + 256 * 256 * r + + end function rgb24bit + + +end module colours_func diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index 2197c3274..4a3414fa6 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -85,7 +85,9 @@ module endfConstants N_da = 117 ,& ! SCONE's fake MT for thermal inelastic scattering N_N_ThermEL = 1002 ,& - N_N_ThermINEL = 1004 + N_N_ThermINEL = 1004 ,& + ! SCONE's fake MT for particle splitting event + N_N_SPLIT = 1005 integer(shortInt),private :: i ! Local, private integer to use array constructor integer(shortInt),parameter :: N_Nl(40) = [(50+i, i =1,40)] diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 9c133c1ee..bf3997446 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -78,7 +78,7 @@ module universalVariables real(defReal), parameter :: joulesPerMeV = 1.60218e-13 ! Convert MeV to J ! Global name variables used to define specific geometry or field types - character(nameLen), parameter :: nameUFS = 'uniFissSites' - character(nameLen), parameter :: nameWW = 'WeightWindows' + character(nameLen), parameter :: nameUFS = 'uniFissSites' + character(nameLen), parameter :: nameWW = 'WeightWindows' end module universalVariables diff --git a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 index 1dbcec26a..85a93fdd6 100644 --- a/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 +++ b/Tallies/TallyClerks/Tests/mgXsClerk_test.f90 @@ -4,8 +4,7 @@ module mgXsClerk_test use endfConstants use genericProcedures, only : numToChar use mgXsClerk_class, only : mgXsClerk - use particle_class, only : particle - use particleDungeon_class, only : particleDungeon + use particle_class, only : particle, particleState use dictionary_class, only : dictionary use scoreMemory_class, only : scoreMemory use testNeutronDatabase_class, only : testNeutronDatabase @@ -104,7 +103,7 @@ subroutine testScoring_clerk1(this) character(:),allocatable :: case type(scoreMemory) :: mem type(particle) :: p - type(particleDungeon) :: pit + type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & nu, chi, P0, P1, P2, P3, P4, & @@ -115,21 +114,24 @@ subroutine testScoring_clerk1(this) call mem % init(1000_longInt, 1) call this % clerk_test1 % setMemAddress(1_longInt) - ! Configure particle dungeon - call pit % init(3) - + ! Report fission particles p % isMG = .false. p % w = 0.5_defReal p % E = 0.3_defReal call p % setMatIdx(2) - call pit % detain(p) + + ! Scoring + pFiss = p + call this % clerk_test1 % reportSpawn(N_FISSION, p, pFiss, this % nucData, mem) p % E = 3.0_defReal - call pit % detain(p) + + ! Scoring + pFiss = p + call this % clerk_test1 % reportSpawn(N_FISSION, p, pFiss, this % nucData, mem) ! Scoring call this % clerk_test1 % reportInColl(p, this % nucData, mem, .false.) - call this % clerk_test1 % reportCycleEnd(pit, mem) p % preCollision % wgt = 0.2_defReal p % preCollision % E = 3.0_defReal @@ -181,7 +183,7 @@ subroutine testScoring_clerk2(this) character(:),allocatable :: case type(scoreMemory) :: mem type(particle) :: p - type(particleDungeon) :: pit + type(particleState) :: pFiss type(outputFile) :: out real(defReal), dimension(:,:), allocatable :: fiss, capt, transFL, transOS, & nu, chi, P0, P1, prod @@ -191,20 +193,24 @@ subroutine testScoring_clerk2(this) call mem % init(1000_longInt, 1) call this % clerk_test2 % setMemAddress(1_longInt) - ! Configure particle dungeon - call pit % init(3) - + ! Report fission particles p % isMG = .false. p % w = 0.5_defReal p % E = 3.0_defReal - call pit % detain(p) + call p % setMatIdx(2) + + ! Scoring + pFiss = p + call this % clerk_test2 % reportSpawn(N_FISSION, p, pFiss, this % nucData, mem) p % E = 0.3_defReal - call pit % detain(p) + + ! Scoring + pFiss = p + call this % clerk_test2 % reportSpawn(N_FISSION, p, pFiss, this % nucData, mem) ! Scoring call this % clerk_test2 % reportInColl(p, this % nucData, mem, .false.) - call this % clerk_test2 % reportCycleEnd(pit, mem) p % preCollision % wgt = 0.2_defReal p % preCollision % E = 0.3_defReal diff --git a/Tallies/TallyClerks/mgXsClerk_class.f90 b/Tallies/TallyClerks/mgXsClerk_class.f90 index 8de338063..fdeda347a 100644 --- a/Tallies/TallyClerks/mgXsClerk_class.f90 +++ b/Tallies/TallyClerks/mgXsClerk_class.f90 @@ -99,7 +99,7 @@ module mgXsClerk_class ! File reports -> run-time procedures procedure :: reportInColl procedure :: reportOutColl - procedure :: reportCycleEnd + procedure :: reportSpawn ! Output procedures procedure :: print @@ -124,7 +124,7 @@ subroutine init(self, dict, name) ! Assign name call self % setName(name) - + ! Load energy map and bin number if (dict % isPresent('energyMap')) then call new_tallyMap(self % energyMap, dict % getDictPtr('energyMap')) @@ -186,7 +186,7 @@ function validReports(self) result(validCodes) class(mgXsClerk),intent(in) :: self integer(shortInt),dimension(:),allocatable :: validCodes - validCodes = [inColl_CODE, outColl_CODE, cycleEnd_CODE] + validCodes = [inColl_CODE, outColl_CODE, spawn_CODE] end function validReports @@ -283,16 +283,16 @@ end subroutine reportInColl !! See tallyClerk_inter for details !! subroutine reportOutColl(self, p, MT, muL, xsData, mem) - class(mgXsClerk), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt), intent(in) :: MT - real(defReal), intent(in) :: muL - class(nuclearDatabase),intent(inout) :: xsData - type(scoreMemory), intent(inout) :: mem - type(particleState) :: preColl, postColl - real(defReal) :: score, prod, mu, mu2, mu3, mu4, mu5 - integer(shortInt) :: enIdx, matIdx, binIdx, binEnOut - integer(longInt) :: addr + class(mgXsClerk), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: MT + real(defReal), intent(in) :: muL + class(nuclearDatabase),intent(inout) :: xsData + type(scoreMemory), intent(inout) :: mem + type(particleState) :: preColl, postColl + real(defReal) :: score, prod, mu, mu2, mu3, mu4, mu5 + integer(shortInt) :: enIdx, matIdx, binIdx, binEnOut + integer(longInt) :: addr ! Get pre and post collision particle state preColl = p % preCollision @@ -406,37 +406,37 @@ subroutine reportOutColl(self, p, MT, muL, xsData, mem) end subroutine reportOutColl !! - !! Process end of the cycle to score fission spectrum with an analog estimator + !! Process fission report !! !! See tallyClerk_inter for details !! - subroutine reportCycleEnd(self, end, mem) - class(mgXsClerk), intent(inout) :: self - class(particleDungeon), intent(in) :: end - type(scoreMemory), intent(inout) :: mem - integer(longInt) :: addr, binIdx, enIdx, matIdx - integer(shortInt) :: N, i + subroutine reportSpawn(self, MT, pOld, pNew, xsData, mem) + class(mgXsClerk), intent(inout) :: self + integer(shortInt), intent(in) :: MT + class(particle), intent(in) :: pOld + class(particleState), intent(in) :: pNew + class(nuclearDatabase), intent(inout) :: xsData + type(scoreMemory), intent(inout) :: mem + integer(longInt) :: addr, binIdx, enIdx, matIdx - ! Loop over the whole neutron population - N = end % popSize() - do i = 1,N + if (MT == N_FISSION) then ! Find bin indexes ! Energy if (allocated(self % energyMap)) then - enIdx = self % energyN + 1 - self % energyMap % map(end % get(i)) + enIdx = self % energyN + 1 - self % energyMap % map(pNew) else enIdx = 1 end if ! Space if (allocated(self % spaceMap)) then - matIdx = self % spaceMap % map(end % get(i)) + matIdx = self % spaceMap % map(pNew) else matIdx = 1 end if ! Return if invalid bin index - if ((enIdx == self % energyN + 1) .or. matIdx == 0) cycle + if ((enIdx == self % energyN + 1) .or. matIdx == 0) return ! Calculate bin address binIdx = self % energyN * (matIdx - 1) + enIdx @@ -445,9 +445,9 @@ subroutine reportCycleEnd(self, end, mem) ! Score energy group of fission neutron call mem % score(ONE, addr + CHI_idx) - end do + end if - end subroutine reportCycleEnd + end subroutine reportSpawn !! !! Final processing to calculate the multi-group cross sections, fission data and diff --git a/Tallies/TallyClerks/tallyClerkSlot_class.f90 b/Tallies/TallyClerks/tallyClerkSlot_class.f90 index bceba3e1a..860ca71b2 100644 --- a/Tallies/TallyClerks/tallyClerkSlot_class.f90 +++ b/Tallies/TallyClerks/tallyClerkSlot_class.f90 @@ -3,7 +3,7 @@ module tallyClerkSlot_class use numPrecision use genericProcedures, only : fatalError use dictionary_class, only : dictionary - use particle_class, only : particle + use particle_class, only : particle, particleState use particleDungeon_class, only : particleDungeon use tallyClerk_inter, only : tallyClerk, setMemAddress_super => setMemAddress, & setName_super => setName, & @@ -45,6 +45,7 @@ module tallyClerkSlot_class procedure :: reportOutColl procedure :: reportPath procedure :: reportTrans + procedure :: reportSpawn procedure :: reportHist procedure :: reportCycleStart procedure :: reportCycleEnd @@ -229,6 +230,24 @@ subroutine reportTrans(self, p, xsData, mem) end subroutine reportTrans + !! + !! Process particle creation report + !! + !! See tallyClerk_inter for details + !! + subroutine reportSpawn(self, MT, pOld, pNew, xsData, mem) + class(tallyClerkSlot), intent(inout) :: self + integer(shortInt), intent(in) :: MT + class(particle), intent(in) :: pOld + class(particleState), intent(in) :: pNew + class(nuclearDatabase), intent(inout) :: xsData + type(scoreMemory), intent(inout) :: mem + + ! Pass call to instance in the slot + call self % slot % reportSpawn(MT, pOld, pNew, xsData, mem) + + end subroutine reportSpawn + !! !! Process history report !! diff --git a/Tallies/TallyClerks/tallyClerk_inter.f90 b/Tallies/TallyClerks/tallyClerk_inter.f90 index 602fea2cb..9eea9eba7 100644 --- a/Tallies/TallyClerks/tallyClerk_inter.f90 +++ b/Tallies/TallyClerks/tallyClerk_inter.f90 @@ -4,7 +4,7 @@ module tallyClerk_inter use tallyCodes use dictionary_class, only : dictionary use genericProcedures, only : fatalError - use particle_class, only : particle + use particle_class, only : particle, particleState use particleDungeon_class, only : particleDungeon use outputFile_class, only : outputFile @@ -57,6 +57,7 @@ module tallyClerk_inter !! reportOutColl -> Process an outgoing from collision report !! reportPath -> Process pathlength report !! reportTrans -> Process transition report + !! reportSpawn -> Process particle generation report !! reportHist -> Process history report !! reportCycleStart -> Process beginning of a cycle report !! reportCycleEnd -> Process end of a cycle report (e.g. Calculate functions of scores like k-eff) @@ -91,6 +92,7 @@ module tallyClerk_inter procedure :: reportOutColl procedure :: reportPath procedure :: reportTrans + procedure :: reportSpawn procedure :: reportHist procedure :: reportCycleStart procedure :: reportCycleEnd @@ -117,6 +119,7 @@ module tallyClerk_inter !! outColl_CODE !! path_CODE !! trans_CODE + !! spawn_CODE !! hist_CODE !! cycleStart_CODE !! cycleEnd_CODE @@ -219,7 +222,7 @@ end subroutine print !! See tallyAdmin_class for implicit assumptions about the report. !! !! Args: - !! p [in] -> Partice + !! p [in] -> Particle !! xsData [inout] -> Nuclear Database with XSs data !! mem [inout] -> Score Memory to put results on !! virtual [in] -> Flag indicating virtual collision @@ -235,7 +238,7 @@ subroutine reportInColl(self,p, xsData, mem, virtual) logical(defBool), intent(in) :: virtual character(100),parameter :: Here = 'reportInColl (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportInColl @@ -246,7 +249,7 @@ end subroutine reportInColl !! See tallyAdmin_class for implicit assumptionas about the report. !! !! Args: - !! p [in] -> Partice + !! p [in] -> Particle !! MT [in] -> MT number of reaction that partilce underwent in the collision !! muL [in] -> Cosine of the collision deflection angle in LAB frame !! xsData [inout]-> Nuclear Database with XSs data @@ -264,7 +267,7 @@ subroutine reportOutColl(self, p, MT, muL, xsData, mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportOutColl (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportOutColl @@ -274,7 +277,7 @@ end subroutine reportOutColl !! See tallyAdmin_class for implicit assumptionas about the report. !! !! Args: - !! p [in] -> Partice + !! p [in] -> Particle !! L [in] -> Length of the path [cm] !! xsData [inout] -> Nuclear Database with XSs data !! mem [inout] -> Score Memory to put results on @@ -290,7 +293,7 @@ subroutine reportPath(self, p, L, xsData,mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportPath (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportPath @@ -300,7 +303,7 @@ end subroutine reportPath !! See tallyAdmin_class for implicit assumptionas about the report. !! !! Args: - !! p [in] -> Partice + !! p [in] -> Particle !! xsData [inout] -> Nuclear Database with XSs data !! mem [inout] -> Score Memory to put results on !! @@ -314,17 +317,45 @@ subroutine reportTrans(self, p, xsData, mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportTrans (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportTrans + !! + !! Process particle creation report + !! + !! See tallyAdmin_class for implicit assumptionas about the report. + !! + !! Args: + !! MT [in] -> MT number of the reaction the particle has undergone + !! pOld [in] -> Particle that caused the branching event + !! pNew [in] -> Particle state of the newly created neutron + !! xsData [inout] -> Nuclear Database with XSs data + !! mem [inout] -> Score Memory to put results on + !! + !! Errors: + !! Depend on specific Clerk + !! + subroutine reportSpawn(self, MT, pOld, pNew, xsData, mem) + class(tallyClerk), intent(inout) :: self + integer(shortInt), intent(in) :: MT + class(particle), intent(in) :: pOld + class(particleState), intent(in) :: pNew + class(nuclearDatabase), intent(inout) :: xsData + type(scoreMemory), intent(inout) :: mem + character(100),parameter :: Here = 'reportSpawn (tallyClerk_inter.f90)' + + call fatalError(Here,'Report was sent to an instance that does not support it.') + + end subroutine reportSpawn + !! !! Process history report !! !! See tallyAdmin_class for implicit assumptionas about the report. !! !! Args: - !! p [in] -> Partice + !! p [in] -> Particle !! xsData [inout] -> Nuclear Database with XSs data !! mem [inout] -> Score Memory to put results on !! @@ -338,7 +369,7 @@ subroutine reportHist(self, p, xsData, mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportHist (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportHist @@ -360,7 +391,7 @@ subroutine reportCycleStart(self, start, mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportCycleStart (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportCycleStart @@ -382,7 +413,7 @@ subroutine reportCycleEnd(self, end, mem) type(scoreMemory), intent(inout) :: mem character(100),parameter :: Here = 'reportCycleEnd (tallyClerk_inter.f90)' - call fatalError(Here,'Report was send to an instance that does not support it.') + call fatalError(Here,'Report was sent to an instance that does not support it.') end subroutine reportCycleEnd diff --git a/Tallies/tallyActiveAdmin_class.f90 b/Tallies/tallyActiveAdmin_class.f90 deleted file mode 100644 index dd7f01d76..000000000 --- a/Tallies/tallyActiveAdmin_class.f90 +++ /dev/null @@ -1,232 +0,0 @@ -module tallyActiveAdmin_class - - use numPrecision - use tallyCodes - use genericProcedures, only : charCmp - use dictionary_class, only : dictionary - use tallyAdminBase_class, only : tallyAdminBase, & - reportInColl_super => reportInColl, & - reportOutColl_super => reportOutColl, & - reportHist_super => reportHist, & - reportCycleStart_super => reportCycleStart, & - reportCycleEnd_super => reportCycleEnd, & - isConverged_super => isConverged, & - init_super => init ,& - print_super => print, & - kill_super => kill ,& - display_super => display - use particle_class, only : particle, phaseCoord - use particleDungeon_class, only : particleDungeon - use keffActiveClerk_class, only : keffActiveClerk - use outputFile_class, only : outputFile - - implicit none - private - - type, public,extends(tallyAdminBase) :: tallyActiveAdmin - private - type(keffActiveClerk) :: keff_estimator - logical(defBool) :: keff_convergence = .false. - contains - ! New Interface - procedure :: keff - - ! Extend superclass procedures - procedure :: reportInColl - procedure :: reportOutColl - procedure :: reportHist - procedure :: reportCycleStart - procedure :: reportCycleEnd - procedure :: init - procedure :: print - procedure :: kill - procedure :: display - procedure :: isConverged - - end type tallyActiveAdmin - -contains - !! - !! Return estimate of k-eff to be used in a calculation for normalisation - !! - function keff(self) result(k) - class(tallyActiveAdmin), intent(in) :: self - real(defReal) :: k - - k = self % keff_estimator % keff() - - end function keff - - !! - !! Process incoming collision report - !! - subroutine reportInColl(self, p) - class(tallyActiveAdmin), intent(inout) :: self - class(particle), intent(in) :: p - - ! Process report with internal Clerk - call self % keff_estimator % reportInColl(p) - - ! Call superclass procedure on self - call reportInColl_super(self, p) - - end subroutine reportInColl - - !! - !! Process outgoing collision report - !! - subroutine reportOutColl(self,p,MT,muL) - class(tallyActiveAdmin), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt), intent(in) :: MT - real(defReal), intent(in) :: muL - - ! Process report with internal Clerk - call self % keff_estimator % reportOutColl(p, MT, muL) - - ! Call superclass procedure on self - call reportOutColl_super(self, p, MT, muL) - - end subroutine reportOutColl - - !! - !! Process history report - !! ASSUMPTIONS: - !! - subroutine reportHist(self, p) - class(tallyActiveAdmin), intent(inout) :: self - class(particle), intent(in) :: p - - ! Process report with internal Clerk - call self % keff_estimator % reportHist(p) - - ! Call superclass procedure on self - call reportHist_super(self, p) - - end subroutine reportHist - - !! - !! Process beginning of a cycle - !! - subroutine reportCycleStart(self, start) - class(tallyActiveAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: start - - ! Process report with internal Clerk - call self % keff_estimator % reportCycleStart(start) - - ! Call superclass procedure on self - call reportCycleStart_super(self, start) - - end subroutine reportCycleStart - - !! - !! Process end of the cycle - !! - subroutine reportCycleEnd(self, end) - class(tallyActiveAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: end - - ! Process report with internal Clerk - call self % keff_estimator % reportCycleEnd(end) - - ! Call superclass procedure on self - call reportCycleEnd_super(self, end) - - end subroutine reportCycleEnd - - !! - !! Initialise active admin - !! - subroutine init(self,dict) - class(tallyActiveAdmin), intent(inout) :: self - class(dictionary),intent(in) :: dict - type(dictionary) :: embDict - character(nameLen) :: entry - real(defReal) :: temp - - call embDict % init(3) - - ! Read settings for embedded clerk and put them into embDict - call dict % getOrDefault(entry,'trigger','no') - - if( charCmp(entry,'yes') ) then - self % keff_convergence = .true. - call dict % get(temp,'SDtarget') - call embDict % store('trigger','yes') - call embDict % store('SDtarget',temp) - end if - - call embDict % store('type','keffActiveClerk') - - ! Initialise embedded clerk - entry = 'A_ADMIN' - call self % keff_estimator % init(embDict,entry) - - ! Load rest of the clerks - call init_super(self,dict) - - end subroutine init - - !! - !! Add all results to outputfile - !! - subroutine print(self,output) - class(tallyActiveAdmin), intent(in) :: self - class(outputFile), intent(inout) :: output - - call self % keff_estimator % print(output) - call print_super(self,output) - - end subroutine print - - !! - !! Deallocates all content - !! - subroutine kill(self) - class(tallyActiveAdmin), intent(inout) :: self - - ! Kill internal Clerks - - ! Call superclass procedure on self - call kill_super(self) - - end subroutine kill - - !! - !! Display results at the end of a cycle - !! - subroutine display(self) - class(tallyActiveAdmin), intent(in) :: self - - print *,'Active Cycle:' - call self % keff_estimator % display() - - ! Call superclass procedure on self - call display_super(self) - - end subroutine display - - !! - !! - !! - function isConverged(self) result(isIt) - class(tallyActiveAdmin), intent(in) :: self - logical(defBool) :: isIt - - if( self % keff_convergence) then - isIt = self % keff_estimator % isConverged() - - if ( self % checkConvergence ) isIt = isIt .and. isConverged_super(self) - - else if ( self % checkConvergence ) then - isIt = isConverged_super(self) - - else - isIt = .false. - - end if - - end function isConverged - -end module tallyActiveAdmin_class diff --git a/Tallies/tallyAdminBase_class.f90 b/Tallies/tallyAdminBase_class.f90 deleted file mode 100644 index 6643e2d4d..000000000 --- a/Tallies/tallyAdminBase_class.f90 +++ /dev/null @@ -1,437 +0,0 @@ -module tallyAdminBase_class - - use numPrecision - use tallyCodes - use genericProcedures, only : fatalError, charCmp - use dictionary_class, only : dictionary - use particle_class, only : particle, phaseCoord - use particleDungeon_class, only : particleDungeon - use tallyClerk_inter, only : tallyClerk - use tallyClerkSlot_class, only : tallyClerkSlot - use tallyClerkFactory_func, only : new_tallyClerk - use outputFile_class, only : outputFile - - - implicit none - private - !! **** MOST LIKLEY CHANGE INTERFACES FOR CLERKS TO INCLUDE FLUX FOR IN COLLISION AND PATH - !! **** PRECALCULATE FLUX HERE SO THERE IS NO NEED TO WARY ABOUT DYNAMIC TYPE OF XSDATA IN CLERKS - !! - !! Base class for the tallies black box. - !! Its responsibilities are as flolow: - !! 1) Accept events reports and routes then to individual tallyClerks - !! 2) Returns k-eff estimate for a current cycle - !! 3) Controls end of calculation - !! 4) Controls printing of calculation progress (to a console) - !! 5) Controls printing of result estimators to a file (filePath and file Format) - !! - !! This class will be extended by inheritance to provide additional functionality - !! i.e. return mesh based weight windows based on fission matrix or similar - !! - !! - type, public:: tallyAdminBase - private - logical(defBool), public :: checkConvergence = .false. - - type(tallyClerkSlot),dimension(:),allocatable :: tallyClerks - - ! Lists of Clerks to be executed for each procedure - integer(shortInt),dimension(:),allocatable :: inCollClerks - integer(shortInt),dimension(:),allocatable :: outCollClerks - integer(shortInt),dimension(:),allocatable :: pathClerks - integer(shortInt),dimension(:),allocatable :: transClerks - integer(shortInt),dimension(:),allocatable :: histClerks - integer(shortInt),dimension(:),allocatable :: cycleStartClerks - integer(shortInt),dimension(:),allocatable :: cycleEndClerks - - ! List of clerks to display - integer(shortInt), dimension(:),allocatable :: displayList - integer(shortInt), dimension(:),allocatable :: triggerList - - contains - ! Report Interface - procedure :: reportInColl - procedure :: reportOutColl - procedure :: reportPath - procedure :: reportTrans - procedure :: reportHist - procedure :: reportCycleStart - procedure :: reportCycleEnd - - ! Display procedures - procedure :: display - - ! Convergance check - procedure :: isConverged - - ! File writing procedures - procedure :: print - - ! Build procedures - procedure :: init - procedure :: addTallyClerk - procedure :: kill - - procedure,private :: addToReports - - end type tallyAdminBase - - - public :: reportInColl - public :: reportOutColl - public :: reportPath - public :: reportTrans - public :: reportHist - public :: reportCycleStart - public :: reportCycleEnd - public :: display - public :: isConverged - public :: init - public :: print - public :: kill - -contains - !! - !! Process incoming collision report - !! - subroutine reportInColl(self,p) - class(tallyAdminBase), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % inCollClerks) -! idx = self % inCollClerks(i) -! call self % tallyClerks(idx) % reportInColl(p) -! -! end do - - end subroutine reportInColl - - !! - !! Display convergance progress of selected tallies on the console - !! - subroutine display(self) - class(tallyAdminBase), intent(in) :: self - integer(shortInt) :: i - -! ! Go through all clerks marked as part of the display -! do i=1,size(self % displayList) -! call self % tallyClerks(i) % display() -! -! end do - - end subroutine display - - !! - !! Perform convergence check in selected clerks - !! - function isConverged(self) result(isIt) - class(tallyAdminBase), intent(in) :: self - logical(defBool) :: isIt - integer(shortInt) :: i,N - -! N = size( self % triggerList) -! -! if( N > 0 ) then -! isIt = self % tallyClerks(1) % isConverged() -! do i = 2,N -! isIt = isIt .and. self % tallyClerks(i) % isConverged() -! -! end do -! else -! isIt = .false. -! -! end if - - end function isConverged - - !! - !! Add all results to outputfile - !! - subroutine print(self,output) - class(tallyAdminBase), intent(in) :: self - class(outputFile), intent(inout) :: output - integer(shortInt) :: i - -! do i=1,size(self % tallyClerks) -! call self % tallyClerks(i) % print(output) -! end do - - end subroutine print - - !! - !! Process outgoing collision report - !! Assume that pre is AFTER any implicit treatment (i.e. implicit capture) - !! - subroutine reportOutColl(self,p,MT,muL) - class(tallyAdminBase), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt), intent(in) :: MT - real(defReal), intent(in) :: muL - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % outCollClerks) -! idx = self % outCollClerks(i) -! call self % tallyClerks(idx) % reportOutColl(p,MT,muL) -! -! end do - - end subroutine reportOutColl - - !! - !! Process pathlength report - !! ASSUMPTIONS: - !! Pathlength must be contained within a single cell and material - !! - subroutine reportPath(self,p,L) - class(tallyAdminBase), intent(inout) :: self - class(particle), intent(in) :: p - real(defReal), intent(in) :: L - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % pathClerks) -! idx = self % pathClerks(i) -! call self % tallyClerks(idx) % reportPath(p,L) -! -! end do - - end subroutine reportPath - - !! - !! Process transition report - !! ASSUMPTIONS: - !! Transition must be a straight line - !! Pre and Post direction is assumed the same (aligned with r_pre -> r_post vector) - !! - subroutine reportTrans(self,p) - class(tallyAdminBase), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % transClerks) -! idx = self % transClerks(i) -! call self % tallyClerks(idx) % reportTrans(p) -! -! end do - - end subroutine reportTrans - - !! - !! Process history report - !! ASSUMPTIONS: - !! - subroutine reportHist(self,p) - class(tallyAdminBase), intent(inout) :: self - class(particle), intent(in) :: p - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % histClerks) -! idx = self % histClerks(i) -! call self % tallyClerks(idx) % reportHist(p) -! -! end do - - - end subroutine reportHist - - !! - !! Process beginning of a cycle - !! - subroutine reportCycleStart(self,start) - class(tallyAdminBase), intent(inout) :: self - class(particleDungeon), intent(in) :: start - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % cycleStartClerks) -! idx = self % cycleStartClerks(i) -! call self % tallyClerks(idx) % reportCycleStart(start) -! -! end do - - end subroutine reportCycleStart - - !! - !! Process end of the cycle - !! - subroutine reportCycleEnd(self,end) - class(tallyAdminBase), intent(inout) :: self - class(particleDungeon), intent(in) :: end - integer(shortInt) :: i, idx - -! ! Go through all clerks that request the report -! do i=1,size(self % cycleEndClerks) -! idx = self % cycleEndClerks(i) -! call self % tallyClerks(idx) % reportCycleEnd(end) -! -! end do - - end subroutine reportCycleEnd - - !! - !! Initialise tallyAdminBase form dictionary - !! - subroutine init(self,dict) - class(tallyAdminBase), intent(inout) :: self - class(dictionary), intent(in) :: dict - character(nameLen),dimension(:),allocatable :: clerks - type(dictionary) :: locDict - character(nameLen) :: entry - logical(defBool) :: partOfDisplay, partOfTriggers - integer(shortInt) :: i - - -! ! Deallocate -! call self % kill() -! -! ! Allocate contents -! allocate(self % inCollClerks(0) ) -! allocate(self % outCollClerks(0) ) -! allocate(self % pathClerks(0) ) -! allocate(self % transClerks(0) ) -! allocate(self % histClerks(0) ) -! allocate(self % cycleStartClerks(0) ) -! allocate(self % cycleEndClerks(0) ) -! -! allocate(self% displayList(0) ) -! allocate(self% triggerList(0) ) -! -! allocate(self % tallyClerks(0)) -! -! ! Read all dictionaries -! call dict % keysDict(clerks) -! -! ! Load all dictionaries as clerks -! do i=1,size(clerks) -! ! Copy dictionary to local copy -! call dict % get(locDict,clerks(i)) -! -! ! Check if it is part of the display -! call locDict % getOrDefault(entry,'display','no') -! partOfDisplay = charCmp(entry,'yes') -! -! ! Check if it is part of the convergance triggers -! call locDict % getOrDefault(entry,'trigger','no') -! partOfTriggers = charCmp(entry,'yes') -! -! ! Get new clerk from factory and store it ina aslot -! call self % addTallyClerk( new_tallyClerk(locDict,clerks(i)),partOfDisplay,partOfTriggers) -! -! end do - - - end subroutine init - - !! - !! Attach new tally - !! - subroutine addTallyClerk(self, clerk, partOfDisplay, partOfTriggers) - class(tallyAdminBase), intent(inout) :: self - class(tallyClerk), intent(in) :: clerk - logical(defBool),intent(in) :: partOfDisplay - logical(defBool),intent(in) :: partOfTriggers - type(tallyClerkSlot) :: localSlot - integer(shortInt),dimension(:),allocatable :: reportCodes - integer(shortInt) :: N, i - - character(100),parameter :: Here = 'addTallyClerk (tallyAdminBase_class.f90)' - -! ! Check if provided clerk is a slot. Give error if it is -! select type(clerk) -! type is (tallyClerkSlot) -! call fatalError(Here,'tallyCleakSlot was passed. It is forbidden to avoid nested slots.') -! end select -! -! ! Check if the tallyAdminBase is initialised -! if( .not. allocated(self % tallyClerks) ) then -! call fatalError(Here,'tallyAdminBase is uninitialised') -! end if -! -! ! Append tally Clerks. Automatic reallocation on assignment. F2008 feature -! localSlot = clerk -! self % tallyClerks = [self % tallyClerks, localSlot] -! -! ! Obtain list of reports requested by the loaded clerk -! N = size(self % tallyClerks) -! reportCodes = self % tallyClerks(N) % validReports() -! -! ! Append report sorting arrays with index of new tallyClerk -! do i=1,size(reportCodes) -! call self % addToReports( reportCodes(i), N ) -! end do -! -! ! If clerk is partOfDisplay append displayList -! if (partOfDisplay) self % displayList = [self % displayList, N ] -! -! ! If clerk is partOfTriggers append triggerList -! if (partOfTriggers) self % triggerList = [self % triggerList, N ] -! if (partOfTriggers) self % checkConvergence = .true. - - end subroutine addTallyClerk - - !! - !! Deallocates all content - !! - subroutine kill(self) - class(tallyAdminBase), intent(inout) :: self - -! if(allocated(self % tallyClerks)) deallocate( self % tallyClerks ) -! -! if(allocated(self % displayList)) deallocate( self % displayList) -! if(allocated(self % triggerList)) deallocate( self % triggerList) -! -! if(allocated(self % inCollClerks)) deallocate( self % inCollClerks) -! if(allocated(self % outCollClerks)) deallocate( self % outCollClerks ) -! if(allocated(self % pathClerks)) deallocate( self % pathClerks ) -! if(allocated(self % transClerks)) deallocate( self % transClerks ) -! if(allocated(self % histClerks)) deallocate( self % histClerks ) -! if(allocated(self % cycleStartClerks)) deallocate( self % cycleStartClerks ) -! if(allocated(self % cycleEndClerks)) deallocate( self % cycleEndClerks ) - - - end subroutine kill - - !! - !! Append sorrting array identified with the code with tallyClerk idx - !! - subroutine addToReports(self,reportCode,idx) - class(tallyAdminBase),intent(inout) :: self - integer(shortInt), intent(in) :: reportCode - integer(shortInt), intent(in) :: idx - character(100),parameter :: Here='addToReports (tallyAdminBase_class.f90)' - -! select case(reportCode) -! case(inColl_CODE) -! self % inCollClerks = [self % inCollClerks, idx] -! -! case(outColl_CODE) -! self % outCollClerks = [ self % outCollClerks, idx] -! -! case(path_CODE) -! self % pathClerks = [ self % pathClerks, idx] -! -! case(trans_CODE) -! self % transClerks = [ self % transClerks, idx] -! -! case(hist_CODE) -! self % histClerks = [ self % histClerks, idx] -! -! case(cycleStart_CODE) -! self % cycleStartClerks = [ self % cycleStartClerks, idx] -! -! case(cycleEnd_CODE) -! self % cycleEndClerks = [ self % cycleEndClerks, idx] -! -! case default -! call fatalError(Here, 'Undefined reportCode') -! end select - - end subroutine addToReports - - -end module tallyAdminBase_class diff --git a/Tallies/tallyAdmin_class.f90 b/Tallies/tallyAdmin_class.f90 index 0e35adfdf..a81b530f0 100644 --- a/Tallies/tallyAdmin_class.f90 +++ b/Tallies/tallyAdmin_class.f90 @@ -6,7 +6,7 @@ module tallyAdmin_class use dictionary_class, only : dictionary use dynArray_class, only : dynIntArray use charMap_class, only : charMap - use particle_class, only : particle + use particle_class, only : particle, particleState use particleDungeon_class, only : particleDungeon use tallyClerk_inter, only : tallyClerk use tallyClerkSlot_class, only : tallyClerkSlot @@ -42,19 +42,20 @@ module tallyAdmin_class !! Private Members: !! atch -> Pointer to an attachment tallyClerk (implements linked-list) !! normBinAddr -> Address of a bin used for normalisation - !! normValue -> Target Value for normalisation - !! normClerkName -> Name of a Clerk used for normalisation - !! tallyClerks -> Array of all defined tally Clerks - !! clerksNameMap -> CharMap that maps Clerk Name to its index in tallyClerks - !! inCollClerks -> List of indices of all Clerks that require inCollReport - !! outCollClerks -> List of indices of all Clerks that require outCollReport - !! pathClerks -> List of indices of all Clerks that require pathReport - !! transClerks -> List of indices of all Clerks that require transReport - !! histClerks -> List of indices of all Clerks that require histReport - !! cycleStartClerks -> List of indices of all Clerks that require cycleStartReport - !! cycleEndClerks -> List of indices of all Clerks that require cycleEndReport - !! displayList -> List of indices of all Clerks that are registered for display - !! mem -> Score Memory for all defined Clerks + !! normValue -> Target value for normalisation + !! normClerkName -> Name of a clerk used for normalisation + !! tallyClerks -> Array of all defined tally clerks + !! clerksNameMap -> CharMap that maps clerk name to its index in tallyClerks + !! inCollClerks -> List of indices of all clerks that require inCollReport + !! outCollClerks -> List of indices of all clerks that require outCollReport + !! pathClerks -> List of indices of all clerks that require pathReport + !! transClerks -> List of indices of all clerks that require transReport + !! spawnClerks -> List of indices of all clerks that require spawnReport + !! histClerks -> List of indices of all clerks that require histReport + !! cycleStartClerks -> List of indices of all clerks that require cycleStartReport + !! cycleEndClerks -> List of indices of all clerks that require cycleEndReport + !! displayList -> List of indices of all clerks that are registered for display + !! mem -> Score Memory for all defined clerks !! !! Interface: !! init -> Initialise from dictionary @@ -66,11 +67,12 @@ module tallyAdmin_class !! reportOutColl -> Process post-collision reports in all clerks !! reportPath -> Process pathlength reports in all clerks !! reportTrans -> Process transition 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 + !! 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 - !! display -> Call "display" on all Clerks registered to display + !! display -> Call "display" on all clerks registered to display !! isConverged -> Return .true. if all convergance targets have been reached !! print -> Prints results to an output file object !! @@ -106,6 +108,7 @@ module tallyAdmin_class type(dynIntArray) :: outCollClerks type(dynIntArray) :: pathClerks type(dynIntArray) :: transClerks + type(dynIntArray) :: spawnClerks type(dynIntArray) :: histClerks type(dynIntArray) :: cycleStartClerks type(dynIntArray) :: cycleEndClerks @@ -131,6 +134,7 @@ module tallyAdmin_class procedure :: reportOutColl procedure :: reportPath procedure :: reportTrans + procedure :: reportSpawn procedure :: reportHist procedure :: reportCycleStart procedure :: reportCycleEnd @@ -263,6 +267,7 @@ recursive subroutine kill(self) call self % outCollClerks % kill() call self % pathClerks % kill() call self % transClerks % kill() + call self % spawnClerks % kill() call self % histClerks % kill() call self % cycleStartClerks % kill() call self % cycleEndClerks % kill() @@ -472,7 +477,7 @@ end subroutine reportInColl !! Process post-collision report !! !! Assumptions: - !! PreCollision state in partice is set to just before this collision + !! PreCollision state in particle is set to just before this collision !! !! Args: !! p [in] -> Particle @@ -583,6 +588,46 @@ recursive subroutine reportTrans(self, p) end subroutine reportTrans + !! + !! Process report for the creation of a new particle + !! + !! Assumptions: + !! It should be sent each time a new particle is created in the simulation + !! by a nuclear reaction or some other mechanism (e.g. splitting) + !! + !! Args: + !! MT [in] -> MT number of the reaction the particle has undergone + !! pOld [in] -> Particle that caused the branching event + !! pNew [in] -> Particle state of the newly created neutron + !! + !! Errors: + !! None + !! + recursive subroutine reportSpawn(self, MT, pOld, pNew) + class(tallyAdmin), intent(inout) :: self + integer(shortInt), intent(in) :: MT + class(particle), intent(in) :: pOld + class(particleState), intent(in) :: pNew + integer(shortInt) :: i, idx + class(nuclearDatabase),pointer :: xsData + character(100), parameter :: Here = "reportSpwan (tallyAdmin_class.f90)" + + ! Call attachment + if(associated(self % atch)) then + call reportSpawn(self % atch, MT, pOld, pNew) + end if + + ! Get Data + xsData => ndReg_get(pOld % getType(), where = Here) + + ! Go through all clerks that request the report + do i=1,self % spawnClerks % getSize() + idx = self % spawnClerks % get(i) + call self % tallyClerks(idx) % reportSpawn(MT, pOld, pNew, xsData, self % mem) + end do + + end subroutine reportSpawn + !! !! Process history report !! @@ -780,6 +825,9 @@ subroutine addToReports(self, reportCode, idx) case(trans_CODE) call self % transClerks % add(idx) + case(spawn_CODE) + call self % spawnClerks % add(idx) + case(hist_CODE) call self % histClerks % add(idx) diff --git a/Tallies/tallyCodes.f90 b/Tallies/tallyCodes.f90 index 44f7bac69..4a05d37aa 100644 --- a/Tallies/tallyCodes.f90 +++ b/Tallies/tallyCodes.f90 @@ -14,9 +14,10 @@ module tallyCodes outColl_CODE = 1001 ,& path_CODE = 1002 ,& trans_CODE = 1003 ,& - hist_CODE = 1004 ,& - cycleStart_CODE = 1005 ,& - cycleEnd_CODE = 1006 + spawn_CODE = 1004 ,& + hist_CODE = 1005 ,& + cycleStart_CODE = 1006 ,& + cycleEnd_CODE = 1007 ! List of codes for fiffrent particle fates integer(shortInt),parameter,public :: abs_FATE = 5000 ,& diff --git a/Tallies/tallyEstimator_class.f90 b/Tallies/tallyEstimator_class.f90 deleted file mode 100644 index c9cf0d36d..000000000 --- a/Tallies/tallyEstimator_class.f90 +++ /dev/null @@ -1,459 +0,0 @@ -module tallyEstimator_class - - use numPrecision - - implicit none - private - - !! - !! Class to contain a score accumulated during a cycle or history or equvalent - !! Internal counter is initialised to 0.0 on creation - !! INTERFACE: - !! add(score) -> adds argument score to internal counter score is defReal,sgortInt or longInt - !! get() -> function, returns value(defReal) of internal counter - !! reset() -> resets internal counter to 0.0 - !! NOTES: - !! This class exists to hide thread safe addition to a tally score from a user. - !! IT IS NOT THREAD SAFE YET. - !! - type, public :: tallyScore - private - real(defReal) :: score = ZERO - contains - generic :: add => add_tallyScore_defReal, & - add_tallyScore_shortInt, & - add_tallyScore_longInt - procedure :: get => get_tallyScore - procedure :: reset => reset_tallyScore - - procedure, private :: add_tallyScore_defReal - procedure, private :: add_tallyScore_shortInt - procedure, private :: add_tallyScore_longInt - - end type tallyScore - - !! - !! Class to store an estimate of a quantity provided over many cycles, histories or batches - !! Is composed and dependent on tallyScore - !! INTERFACE: - !! addEstimate(est) -> adds extra estimate to internal storage, estimate(est) can be - !! defReal, shortInt or longInt. - !! getEstimate(est,STD,N)-> calculates mean estimate and standard deviation of an estimate. - !! est & STD are defReal. N is shortInt. There is no protection - !! against N = 0. Compiler default rules for division by 0 are - !! followed. For N=1 pseudo-variance estimate is returned with N-1 - !! term dropped. This is for estetic reasons. This estimate is - !! unreliable. - !! reset() -> resets internal counters to 0.0 - !! NOTES: - !! This class exists to wrap cumulative sum and of estimate and estimate^2. - !! It is assumed it will be called beween parallel cycles. IS NOT THREAD SAFE - !! - type, public :: tallyCounter - private - type(tallyScore) :: cSum - type(tallyScore) :: cSum2 - - contains - generic :: addEstimate => addEstimate_tallyCounter_defReal, & - addEstimate_tallyCounter_shortInt, & - addEstimate_tallyCounter_longInt - generic :: getEstimate => getEstimate_tallyCounter_withSTD, & - getEstimate_tallyCounter_withoutSTD - procedure :: reset => reset_tallyCounter - - procedure,private :: getEstimate_tallyCounter_withSTD - procedure,private :: getEstimate_tallyCounter_withoutSTD - - procedure,private :: addEstimate_tallyCounter_defReal - procedure,private :: addEstimate_tallyCounter_shortInt - procedure,private :: addEstimate_tallyCounter_longInt - - end type tallyCounter - - !! - !! Class that combine tallyScore and tally Counter into one entity. - !! It is used with estimates based on a single integration (i.e. reaction rates) - !! INTERFACE: - !! add(score) -> adds argument score to internal tallyScore is defReal,sgortInt or longInt - !! get() -> function, returns value(defReal) in internal tallyScore - !! closeBatch(N)-> closes batch. Normalises value in internal tallyScore by N. - !! N can be defReal, shortInt and longInt. If no scores were accumulated - !! during a batch (exactly 0.0 value in tallyScore), batch-wise estimate is 0.0 - !! even if N=0. Otherwise Compiler default rules for division by 0 are followed. - !! getEstimate(est,STD,Nb)-> calculates mean estimate and standard deviation of an estimate. - !! Assumes that Nb batches vere closed before. - !! est & STD are defReal. Nb is shortInt. There is no protection - !! against Nb = 0. Compiler default rules for division by 0 are - !! followed. For Nb=1 pseudo-variance estimate is returned with Nb-1 - !! term dropped. This is for estetic reasons. This estimate is - !! unreliable. - !! NOTES: - !! Class does not store how many batches were closed to save memory and avoid redundunt - !! additions in each class instance. - !! - type, public :: tallyEstimator - private - type(tallyScore) :: score - type(tallyCounter) :: estimate - contains - generic :: add => add_tallyEstimator_defReal, & - add_tallyEstimator_shortInt, & - add_tallyEstimator_longInt - procedure :: get => get_tallyEstimator - generic :: closeBatch => closeBatch_tallyEstimator_defReal, & - closeBatch_tallyEstimator_shortInt, & - closeBatch_tallyEstimator_longInt - generic :: getEstimate => getEstimate_tallyEstimator_withSTD, & - getEstimate_tallyEstimator_withoutSTD - procedure :: reset => reset_tallyEstimator - - procedure, private :: add_tallyEstimator_defReal - procedure, private :: add_tallyEstimator_shortInt - procedure, private :: add_tallyEstimator_longInt - - procedure, private :: closeBatch_tallyEstimator_defReal - procedure, private :: closeBatch_tallyEstimator_shortInt - procedure, private :: closeBatch_tallyEstimator_longInt - - procedure, private :: getEstimate_tallyEstimator_withSTD - procedure, private :: getEstimate_tallyEstimator_withoutSTD - - end type tallyEstimator -contains - -!**************************************************************************************************! -!* tallyScore class procedures *! -!**************************************************************************************************! - !! - !! Add score(defReal) to a score counter - !! - elemental subroutine add_tallyScore_defReal(self,score) - class(tallyScore), intent(inout) :: self - real(defReal), intent(in) :: score - - self % score = self % score + score - - end subroutine add_tallyScore_defReal - - !! - !! Add score(shortInt) to a score counter - !! - elemental subroutine add_tallyScore_shortInt(self,score) - class(tallyScore), intent(inout) :: self - integer(shortInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert to real - score_tmp = real(score,defReal) - - ! Call procedure for real - call self % add_tallyScore_defReal(score_tmp) - - end subroutine add_tallyScore_shortInt - - !! - !! Add score(longInt) to a score counter - !! - elemental subroutine add_tallyScore_longInt(self,score) - class(tallyScore), intent(inout) :: self - integer(longInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert to real - score_tmp = real(score,defReal) - - ! Call procedure for real - call self % add_tallyScore_defReal(score_tmp) - - end subroutine add_tallyScore_longInt - - !! - !! Returns accumulated score value - !! - elemental function get_tallyScore(self) result(score) - class(tallyScore), intent(in) :: self - real(defReal) :: score - - score = self % score - - end function get_tallyScore - - !! - !! Resets score counter - !! - elemental subroutine reset_tallyScore(self) - class(tallyScore), intent(inout) :: self - - self % score = ZERO - - end subroutine reset_tallyScore - -!**************************************************************************************************! -!* tallyCounter class procedures *! -!**************************************************************************************************! - !! - !! Add a result estimate(defReal) to a tallyCounter - !! - elemental subroutine addEstimate_tallyCounter_defReal(self,score) - class(tallyCounter), intent(inout) :: self - real(defReal),intent(in) :: score - - call self % cSum % add(score) - call self % cSum2 % add(score * score ) - - end subroutine addEstimate_tallyCounter_defReal - - !! - !! Add a result estimate(shortInt) to a tallyCounter - !! - elemental subroutine addEstimate_tallyCounter_shortInt(self,score) - class(tallyCounter), intent(inout) :: self - integer(shortInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert and call procedure for defReal - score_tmp = real(score,defReal) - call self % addEstimate(score_tmp) - - end subroutine addEstimate_tallyCounter_shortInt - - !! - !! Add a result estimate(longInt) to a tallyCounter - !! - elemental subroutine addEstimate_tallyCounter_longInt(self,score) - class(tallyCounter), intent(inout) :: self - integer(longInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert and call procedure for defReal - score_tmp = real(score,defReal) - call self % addEstimate(score_tmp) - - end subroutine addEstimate_tallyCounter_longInt - - !! - !! Obtain estimate from a tally Counter. Returns standard deviation - !! Normalise assuming N (shortInt) samples were collected - !! There is no protection against -ve and N = 0 normal float arthmetic rules are followed - !! - elemental subroutine getEstimate_tallyCounter_withSTD(self,est,STD,N) - class(tallyCounter), intent(in) :: self - real(defReal), intent(out) :: est - real(defReal), intent(out) :: STD - integer(shortInt), intent(in) :: N - real(defReal) :: cSum2 - integer(shortInt) :: Nm1 - - cSum2 = self % cSum2 % get() - - ! Call getEstimate without STD to calculate mean - call self % getEstimate(est,N) - - ! Protect against STD = NaN if N == 1. Return variance equal to estimate - ! This mostly done for aesthetics. Single sample estimate is unrealible - if ( N == 1) then - STD = est - - else ! when N /= 1 - ! Precalculate denominator (N-1) - Nm1 = N - 1 - - ! Calculate the Variance of the mean - STD = cSum2 / Nm1 /N - est * est / Nm1 - - ! Calculate Standard Deviation - STD = sqrt(STD) - - end if - - end subroutine getEstimate_tallyCounter_withSTD - - !! - !! Obtain estimate from a tally Counter. Does not returns Standard Deviation. - !! Normalise assuming N (shortInt) samples were collected - !! There is no protection against -ve and N = 0 normal float arthmetic rules are followed - !! - elemental subroutine getEstimate_tallyCounter_withoutSTD(self,est,N) - class(tallyCounter), intent(in) :: self - real(defReal), intent(out) :: est - integer(shortInt), intent(in) :: N - real(defReal) :: cSum - - cSum = self % cSum % get() - - ! Calculate Mean - est = cSum / N - - end subroutine getEstimate_tallyCounter_withoutSTD - - !! - !! Reset tallyCounter - !! - elemental subroutine reset_tallyCounter(self) - class(tallyCounter), intent(inout) :: self - - call self % cSum % reset() - call self % cSum2 % reset() - - end subroutine reset_tallyCounter - -!**************************************************************************************************! -!* tallyEstimate class procedures *! -!**************************************************************************************************! - - !! - !! Add new score(defReal) to a tallyScore in the estimator - !! - elemental subroutine add_tallyEstimator_defReal(self,score) - class(tallyEstimator), intent(inout) :: self - real(defReal), intent(in) :: score - - call self % score % add(score) - - end subroutine add_tallyEstimator_defReal - - !! - !! Add new score(defReal) to a tallyScore in the estimator - !! - elemental subroutine add_tallyEstimator_shortInt(self,score) - class(tallyEstimator), intent(inout) :: self - integer(shortInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert and call procedure for defReal - score_tmp = real(score,defReal) - call self % score % add(score) - - end subroutine add_tallyEstimator_shortInt - - !! - !! Add new score(defReal) to a tallyScore in the estimator - !! - elemental subroutine add_tallyEstimator_longInt(self,score) - class(tallyEstimator), intent(inout) :: self - integer(longInt), intent(in) :: score - real(defReal) :: score_tmp - - ! Convert and call procedure for defReal - score_tmp = real(score,defReal) - call self % score % add(score) - - end subroutine add_tallyEstimator_longInt - - !! - !! Return accumulated score - !! - elemental function get_tallyEstimator(self) result (score) - class(tallyEstimator), intent(in) :: self - real(defReal) :: score - - score = self % score % get() - - end function get_tallyEstimator - - !! - !! Closes Batch - !! Normalise score assuming that total weight of the samples was N(defReal) - !! There is no protection against N == 0 - !! normal float arthmetic rules are followed - !! Special case: if score == 0.0 (exactly in floating points) - !! It corresponds (most likley) to no samples beeing collected - !! Then batch-wise mean is estimated to be 0.0 - !! - elemental subroutine closeBatch_tallyEstimator_defReal(self,N) - class(tallyEstimator), intent(inout) :: self - real(defReal), intent(in) :: N - real(defReal) :: score, mean - - ! Extract score and reset tallyScore - score = self % score % get() - call self % score % reset() - - ! Detect spetial case when no samples were collected - if( score == ZERO ) then - mean = ZERO - - else - mean = score / N - - end if - - ! Add estimate of the mean to the tallyCounter - call self % estimate % addEstimate(mean) - - end subroutine - - !! - !! Normalise score assuming that total weight of the samples was N(shortInt) - !! Converts N into defReal and calls procedure for defReal - !! - elemental subroutine closeBatch_tallyEstimator_shortInt(self,N) - class(tallyEstimator), intent(inout) :: self - integer(shortInt), intent(in) :: N - real(defReal) :: N_tmp - - ! Convert and call procedure for defReal - N_tmp = real(N,defReal) - call self % closeBatch(N_tmp) - - end subroutine closeBatch_tallyEstimator_shortInt - - !! - !! Normalise score assuming that total weight of the samples was N(longInt) - !! Converts N into defReal and calls procedure for defReal - !! - elemental subroutine closeBatch_tallyEstimator_longInt(self,N) - class(tallyEstimator), intent(inout) :: self - integer(longInt), intent(in) :: N - real(defReal) :: N_tmp - - ! Convert and call procedure for defReal - N_tmp = real(N,defReal) - call self % closeBatch(N_tmp) - - end subroutine closeBatch_tallyEstimator_longInt - - !! - !! Obtain estimate from a tally Counter. Returns Standard Deviation. - !! Normalise assuming N (shortInt) batches were collected (with closeBatch) - !! There is no protection against N < 0 and N == 0 - !! Normal float arthmetic rules are followed - !! - elemental subroutine getEstimate_tallyEstimator_withSTD(self,est,STD,Nb) - class(tallyEstimator), intent(in) :: self - real(defReal), intent(out) :: est - real(defReal), intent(out) :: STD - integer(shortInt), intent(in) :: Nb - - call self % estimate % getEstimate(est,STD,Nb) - - end subroutine getEstimate_tallyEstimator_withSTD - - !! - !! Obtain estimate from a tally Counter. Does not return Standard Deviation. - !! Normalise assuming N (shortInt) batches were collected (with closeBatch) - !! There is no protection against N < 0 and N == 0 - !! Normal float arthmetic rules are followed - !! - elemental subroutine getEstimate_tallyEstimator_withoutSTD(self,est,Nb) - class(tallyEstimator), intent(in) :: self - real(defReal), intent(out) :: est - integer(shortInt), intent(in) :: Nb - - call self % estimate % getEstimate(est,Nb) - - end subroutine getEstimate_tallyEstimator_withoutSTD - - !! - !! Resets score and estimate of tallyEstimate - !! - elemental subroutine reset_tallyEstimator(self) - class(tallyEstimator), intent(inout) :: self - - call self % score % reset() - call self % estimate % reset() - - end subroutine reset_tallyEstimator - -end module tallyEstimator_class diff --git a/Tallies/tallyInactiveAdmin_class.f90 b/Tallies/tallyInactiveAdmin_class.f90 deleted file mode 100644 index 1f7e4507a..000000000 --- a/Tallies/tallyInactiveAdmin_class.f90 +++ /dev/null @@ -1,132 +0,0 @@ -module tallyInactiveAdmin_class - - use numPrecision - use tallyCodes - use dictionary_class, only : dictionary - use tallyAdminBase_class, only : tallyAdminBase, & - reportCycleStart_super => reportCycleStart, & - reportCycleEnd_super => reportCycleEnd, & - init_super => init, & - kill_super => kill, & - display_super => display - use particle_class, only : particle, phaseCoord - use particleDungeon_class, only : particleDungeon - use keffInactiveClerk_class, only : keffInactiveClerk - - implicit none - private - - type, public,extends(tallyAdminBase) :: tallyInactiveAdmin - private - type(keffInactiveClerk) :: keff_estimator - contains - ! New interface - procedure :: keff - - ! Extend superclass procedures - procedure :: reportCycleStart - procedure :: reportCycleEnd - procedure :: init - procedure :: kill - procedure :: display - - ! Overwrite superclass procedures - - - end type tallyInactiveAdmin - -contains - - !! - !! Return estimate of k-eff to be used in a calculation for normalisation - !! - function keff(self) result(k) - class(tallyInactiveAdmin), intent(in) :: self - real(defReal) :: k - - k = self % keff_estimator % keff() - - end function keff - - !! - !! Process beginning of a cycle - !! - subroutine reportCycleStart(self,start) - class(tallyInactiveAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: start - - ! Process report with internal Clerk - call self % keff_estimator % reportCycleStart(start) - - ! Call superclass procedure on self - call reportCycleStart_super(self,start) - - end subroutine reportCycleStart - - !! - !! Process end of the cycle - !! - subroutine reportCycleEnd(self,end) - class(tallyInactiveAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: end - - ! Process report with internal Clerk - call self % keff_estimator % reportCycleEnd(end) - - ! Call superclass procedure on self - call reportCycleEnd_super(self,end) - - end subroutine reportCycleEnd - - !! - !! Initialise active admin - !! - subroutine init(self,dict) - class(tallyInactiveAdmin), intent(inout) :: self - class(dictionary),intent(in) :: dict - type(dictionary) :: embDict - character(nameLen) :: name - - ! Get settings for the embedded clerk into dictionary - call embDict % init(1) - call embDict % store('type','keffInactiveClerk') - - ! Initialise embedded clerk - name = 'IA_ADMIN' - call self % keff_estimator % init(embDict,name) - - ! Load rest of the clerks - call init_super(self,dict) - - end subroutine init - - - - !! - !! Deallocates all content - !! - subroutine kill(self) - class(tallyInactiveAdmin), intent(inout) :: self - - ! Kill internal Clerks - - ! Call superclass procedure on self - call kill_super(self) - - end subroutine kill - - !! - !! Display results at the end of a cycle - !! - subroutine display(self) - class(tallyInactiveAdmin), intent(in) :: self - - print *,'Inactive Cycle:' - call self % keff_estimator % display() - - ! Call superclass procedure on self - call display_super(self) - - end subroutine display - -end module tallyInactiveAdmin_class diff --git a/Tallies/tallyTimeAdmin_class.f90 b/Tallies/tallyTimeAdmin_class.f90 deleted file mode 100644 index 25163d72f..000000000 --- a/Tallies/tallyTimeAdmin_class.f90 +++ /dev/null @@ -1,219 +0,0 @@ -module tallyTimeAdmin_class - - use numPrecision - use tallyCodes - use genericProcedures, only : charCmp - use dictionary_class, only : dictionary - use tallyAdminBase_class, only : tallyAdminBase, & - reportInColl_super => reportInColl, & - reportHist_super => reportHist, & - reportCycleStart_super => reportCycleStart, & - reportCycleEnd_super => reportCycleEnd, & - isConverged_super => isConverged, & - init_super => init ,& - print_super => print, & - kill_super => kill ,& - display_super => display - use particle_class, only : particle, phaseCoord - use particleDungeon_class, only : particleDungeon - use timeClerk_class, only : timeClerk - use outputFile_class, only : outputFile - - implicit none - private - - type, public,extends(tallyAdminBase) :: tallyTimeAdmin - private - type(timeClerk) :: power_estimator - logical(defBool):: power_convergence = .false. - contains - ! New Interface - procedure :: power - procedure :: stepLength - procedure :: incrementStep - - ! Extend superclass procedures - procedure :: reportInColl - procedure :: reportHist - procedure :: reportCycleStart - procedure :: reportCycleEnd - procedure :: init - procedure :: print - procedure :: kill - procedure :: display - procedure :: isConverged - - end type tallyTimeAdmin - -contains - !! - !! Return estimate of power to be used in normalisation and results - !! - function power(self) result(p) - class(tallyTimeAdmin), intent(in) :: self - real(defReal) :: p - - p = self % power_estimator % power() - - end function power - - !! - !! Return timestep length given an index - !! - function stepLength(self,idx) result(dt) - class(tallyTimeAdmin), intent(in) :: self - integer(shortInt), intent(in) :: idx - real(defReal) :: dt - - dt = self % power_estimator % stepLength(idx) - - end function stepLength - - !! - !! Increment step count - !! - subroutine incrementStep(self) - class(tallyTimeAdmin), intent(inout) :: self - call self % power_estimator % incrementStep() - end subroutine incrementStep - - !! - !! Process incoming collision report - !! - subroutine reportInColl(self, p) - class(tallyTimeAdmin), intent(inout) :: self - class(particle), intent(in) :: p - - ! Process report with internal Clerk - call self % power_estimator % reportInColl(p) - - ! Call superclass procedure on self - call reportInColl_super(self, p) - - end subroutine reportInColl - - !! - !! Process history report - !! ASSUMPTIONS: - !! - subroutine reportHist(self, p) - class(tallyTimeAdmin), intent(inout) :: self - class(particle), intent(in) :: p - - ! Process report with internal Clerk - call self % power_estimator % reportHist(p) - - ! Call superclass procedure on self - call reportHist_super(self, p) - - end subroutine reportHist - - !! - !! Process beginning of a cycle - !! - subroutine reportCycleStart(self, start) - class(tallyTimeAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: start - - ! Process report with internal Clerk - call self % power_estimator % reportCycleStart(start) - - ! Call superclass procedure on self - call reportCycleStart_super(self, start) - - end subroutine reportCycleStart - - !! - !! Process end of the cycle - !! - subroutine reportCycleEnd(self, end) - class(tallyTimeAdmin), intent(inout) :: self - class(particleDungeon), intent(in) :: end - - ! Process report with internal Clerk - call self % power_estimator % reportCycleEnd(end) - - ! Call superclass procedure on self - call reportCycleEnd_super(self, end) - - end subroutine reportCycleEnd - - !! - !! Add all results to outputfile - !! - subroutine print(self,output) - class(tallyTimeAdmin), intent(in) :: self - class(outputFile), intent(inout) :: output - - print *, 'HERE' - call self % power_estimator % print(output) - call print_super(self,output) - - end subroutine print - - !! - !! Initialise time admin - !! - subroutine init(self,dict) - class(tallyTimeAdmin), intent(inout) :: self - class(dictionary),intent(in) :: dict - character(nameLen) :: entry - - ! Initialise embedded clerk - entry = 'T_ADMIN' - call self % power_estimator % init(dict,entry) - - ! Load rest of the clerks - call init_super(self,dict) - - end subroutine init - - !! - !! Deallocates all content - !! - subroutine kill(self) - class(tallyTimeAdmin), intent(inout) :: self - - ! Kill internal Clerks - - ! Call superclass procedure on self - call kill_super(self) - - end subroutine kill - - !! - !! Display results at the end of a cycle - !! - subroutine display(self) - class(tallyTimeAdmin), intent(in) :: self - - call self % power_estimator % display() - - ! Call superclass procedure on self - call display_super(self) - - end subroutine display - - !! - !! - !! - function isConverged(self) result(isIt) - class(tallyTimeAdmin), intent(in) :: self - logical(defBool) :: isIt - - if( self % power_convergence) then - isIt = self % power_estimator % isConverged() - - if ( self % checkConvergence ) isIt = isIt .and. isConverged_super(self) - - else if ( self % checkConvergence ) then - isIt = isConverged_super(self) - - else - isIt = .false. - - end if - - end function isConverged - -end module tallyTimeAdmin_class diff --git a/UserInterface/fileOutput/CMakeLists.txt b/UserInterface/fileOutput/CMakeLists.txt index bf913adf4..3ad697235 100644 --- a/UserInterface/fileOutput/CMakeLists.txt +++ b/UserInterface/fileOutput/CMakeLists.txt @@ -4,6 +4,7 @@ add_sources( ./outputFile_class.f90 ./asciiOutputFactory_func.f90 ./asciiMATLAB_class.f90 ./asciiJSON_class.f90 - ./dummyPrinter_class.f90) + ./dummyPrinter_class.f90 + ./delayedStream_class.f90) add_unit_tests (./Tests/outputFile_test.f90) diff --git a/UserInterface/fileOutput/asciiJSON_class.f90 b/UserInterface/fileOutput/asciiJSON_class.f90 index 18ddf36eb..4fd6fa66c 100644 --- a/UserInterface/fileOutput/asciiJSON_class.f90 +++ b/UserInterface/fileOutput/asciiJSON_class.f90 @@ -21,7 +21,7 @@ module asciiJSON_class !! There is indentation for each block and entry, but long multi-dimensional (N-D) arrays !! are printed on a single line. Also N-D arrays are represented as a JSON array of arrays. !! - !! JSON output is fairly standard so it can be easlily read into Python and other enviroments + !! JSON output is fairly standard so it can be easily read into Python and other environments !! !! TODO: !! Currently there is a dirty fix to remove commas at the end of a block, which requires to @@ -39,7 +39,6 @@ module asciiJSON_class !! type, public, extends(asciiOutput) :: asciiJSON private - type(charTape) :: output integer(shortInt) :: ind_lvl = 0 integer(shortInt), dimension(:), allocatable :: shapeBuffer @@ -49,8 +48,8 @@ module asciiJSON_class contains procedure :: init + procedure :: endFile procedure :: extension - procedure :: writeToFile procedure :: startBlock procedure :: endBlock @@ -71,51 +70,43 @@ subroutine init(self) class(asciiJSON), intent(inout) :: self ! Add the initial bracket - call self % output % append( "{" // NEWLINE) + call self % append( "{" // NEWLINE) self % ind_lvl = 1 end subroutine init !! - !! Return approperiate extension for the file + !! Finalise the printer !! !! See asciiOutput_inter for details !! - pure function extension(self) result(str) - class(asciiJSON), intent(in) :: self - character(:), allocatable :: str - - str = 'json' - - end function extension - - !! - !! Print the output to the given unit - !! - !! See asciiOutput_inter for details - !! - subroutine writeToFile(self, unit) + subroutine endFile(self) class(asciiJSON), intent(inout) :: self - integer(shortInt), intent(in) :: unit - character(:),allocatable :: form ! Dirty fix ! Remove comma from the last entry by rewind ! Need to check that the character is a comma to avoid removing { when there is an empty block ! TODO: Find better way. Will clash with any proper stream output - if (self % output % get(self % output % length() - 1) == ",") then - call self % output % cut(2) - call self % output % append(NEWLINE) + if (self % peek(2) == ",") then + call self % cut(2) + call self % append(NEWLINE) end if + call self % append("}" // NEWLINE) - ! Close the file - call self % output % append("}") + end subroutine endFile - form = '(A' // numToChar(self % output % length()) // ')' + !! + !! Return appropriate extension for the file + !! + !! See asciiOutput_inter for details + !! + pure function extension(self) result(str) + class(asciiJSON), intent(in) :: self + character(:), allocatable :: str - write(unit,form) self % output % expose() + str = 'json' - end subroutine writeToFile + end function extension !! !! Change state to writing new block with "name" @@ -127,10 +118,10 @@ subroutine startBlock(self, name) character(nameLen), intent(in) :: name ! Write indentation - call self % output % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) + call self % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) ! Open new block and increase indentation - call self % output % append(QUOTE // trim(name) // QUOTE // ":" // "{" // NEWLINE) + call self % append(QUOTE // trim(name) // QUOTE // ":" // "{" // NEWLINE) self % ind_lvl = self % ind_lvl + 1 end subroutine startBlock @@ -147,17 +138,17 @@ subroutine endBlock(self) ! Remove comma from the last entry by rewind ! Need to check that the character is a comma to avoid removing { when there is an empty block ! TODO: Find better way. Will clash with any proper stream output - if (self % output % get(self % output % length() - 1) == ",") then - call self % output % cut(2) - call self % output % append(NEWLINE) + if (self % peek(2) == ",") then + call self % cut(2) + call self % append(NEWLINE) end if ! Decrease and write indentation self % ind_lvl = self % ind_lvl - 1 - call self % output % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) + call self % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) ! Close the block - call self % output % append("}" // ","// NEWLINE) + call self % append("}" // ","// NEWLINE) end subroutine endBlock @@ -171,10 +162,10 @@ subroutine startEntry(self, name) character(*), intent(in) :: name ! Print indentation - call self % output % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) + call self % append(repeat(BLANK, self % ind_lvl * IND_SIZE)) ! Write the name - call self % output % append(QUOTE // trim(name) // QUOTE // ":") + call self % append(QUOTE // trim(name) // QUOTE // ":") end subroutine startEntry @@ -188,7 +179,7 @@ subroutine endEntry(self) ! End with NEWLINE ! Comma will be written together with a number/char - call self % output % append(NEWLINE) + call self % append(NEWLINE) end subroutine endEntry @@ -235,12 +226,12 @@ subroutine printNum(self, val) mul = 1 do i = 1, size(self % shapeBuffer) mul = mul * self % shapeBuffer(i) - if (modulo(self % count, mul) == 0) call self % output % append("[") + if (modulo(self % count, mul) == 0) call self % append("[") end do end if ! Print the number - call self % output % append(val) + call self % append(val) if (self % in_array) then ! Append the count @@ -250,12 +241,12 @@ subroutine printNum(self, val) mul = 1 do i = 1, size(self % shapeBuffer) mul = mul * self % shapeBuffer(i) - if (modulo(self % count, mul) == 0) call self % output % append("]") + if (modulo(self % count, mul) == 0) call self % append("]") end do end if ! Finish entry with a comma - call self % output % append(",") + call self % append(",") end subroutine printNum @@ -274,12 +265,12 @@ subroutine printChar(self, val) mul = 1 do i = 1, size(self % shapeBuffer) mul = mul * self % shapeBuffer(i) - if (modulo(self % count, mul) == 0) call self % output % append("[") + if (modulo(self % count, mul) == 0) call self % append("[") end do end if ! Print the character - call self % output % append(QUOTE // val // QUOTE) + call self % append(QUOTE // val // QUOTE) if (self % in_array) then ! Append the count @@ -289,12 +280,12 @@ subroutine printChar(self, val) mul = 1 do i = 1, size(self % shapeBuffer) mul = mul * self % shapeBuffer(i) - if (modulo(self % count, mul) == 0) call self % output % append("]") + if (modulo(self % count, mul) == 0) call self % append("]") end do end if ! Finish entry with a comma - call self % output % append(",") + call self % append(",") end subroutine printChar diff --git a/UserInterface/fileOutput/asciiMATLAB_class.f90 b/UserInterface/fileOutput/asciiMATLAB_class.f90 index eeb232789..b844dc710 100644 --- a/UserInterface/fileOutput/asciiMATLAB_class.f90 +++ b/UserInterface/fileOutput/asciiMATLAB_class.f90 @@ -23,8 +23,8 @@ module asciiMATLAB_class !! !! Printer for ASCII MATLAB output file !! - !! Creates a computer-readable matalab ".m" file. - !! Despite beeing in ASCII it is not intended for human-readability: + !! Creates a computer-readable matlab ".m" file. + !! Despite being in ASCII it is not intended for human-readability: !! -> Each entry is a single line. !! -> Furthermore every array is printed as RANK 1 array and a reshape function. !! @@ -40,7 +40,7 @@ module asciiMATLAB_class !! blockLevel -> Current block level !! output -> Character that contain the output !! prefix -> Current prefix for the block - !! shapeBuffer -> Buffored shape of the current array + !! shapeBuffer -> Buffered shape of the current array !! !! Interface: !! asciiOutput Interface @@ -52,17 +52,14 @@ module asciiMATLAB_class type(stackChar) :: blockNameStack integer(shortInt) :: blockLevel = 0 - ! Outputfile - type(charTape) :: output - ! Buffers type(charTape) :: prefix integer(shortInt), dimension(:), allocatable :: shapeBuffer contains procedure :: init + procedure :: endFile procedure :: extension - procedure :: writeToFile procedure :: startBlock procedure :: endBlock @@ -90,33 +87,29 @@ subroutine init(self) end subroutine init !! - !! Return approperiate extension for the file + !! Finalise the printer !! !! See asciiOutput_inter for details !! - pure function extension(self) result(str) - class(asciiMATLAB), intent(in) :: self - character(:), allocatable :: str + subroutine endFile(self) + class(asciiMATLAB), intent(inout) :: self - str = 'm' + ! Nothing to do - end function extension + end subroutine endFile !! - !! Print the output to the given unit + !! Return appropriate extension for the file !! !! See asciiOutput_inter for details !! - subroutine writeToFile(self, unit) - class(asciiMATLAB), intent(inout) :: self - integer(shortInt), intent(in) :: unit - character(:),allocatable :: form - - form = '(A' // numToChar(self % output % length()) // ')' + pure function extension(self) result(str) + class(asciiMATLAB), intent(in) :: self + character(:), allocatable :: str - write(unit,form) self % output % expose() + str = 'm' - end subroutine writeToFile + end function extension !! !! Change state to writing new block with "name" @@ -128,11 +121,6 @@ subroutine startBlock(self, name) character(nameLen), intent(in) :: name character(100), parameter :: Here ='startBlock (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state == IN_ENTRY .or. self % state == IN_ARRAY) then - ! call fatalError(Here,'Cannot start writing new block inside entry or array') - ! end if - ! Update state - change current prefix call self % blockNameStack % push(name) call self % prefix % append(trim(name) // '_') @@ -151,15 +139,6 @@ subroutine endBlock(self) integer(shortInt) :: N character(100), parameter :: Here ='endBlock (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state == IN_ENTRY .or. self % state == IN_ARRAY) then - ! call fatalError(Here,'Cannot end writing new block inside entry or array') - ! end if - ! - ! if ( self % blockLevel == 0) then - ! call fatalError(Here,'Cannot exit from root block') - ! end if - ! Update state - change current prefix call self % blockNameStack % pop(temp) @@ -180,17 +159,12 @@ subroutine startEntry(self, name) character(*), intent(in) :: name character(100), parameter :: Here ='startEntry (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state == IN_ENTRY .or. self % state == IN_ARRAY) then - ! call fatalError(Here,'Cannot star writing new entry inside entry or array') - ! end if - ! Update state self % state = IN_ENTRY ! Write variable name with prefix - call self % output % append( trim(self % prefix % expose()) // trim(name)) - call self % output % append( ' = ') + call self % append( trim(self % prefix % expose()) // trim(name)) + call self % append( ' = ') end subroutine startEntry @@ -204,16 +178,11 @@ subroutine endEntry(self) class(asciiMATLAB), intent(inout) :: self character(100), parameter :: Here ='endEntry (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state == IN_BLOCK .or. self % state == IN_ARRAY) then - ! call fatalError(Here,'Cannot finish entry in block or array') - ! end if - ! Update state self % state = IN_BLOCK ! Write semicolon and newline - call self % output % append( ';' // NEWLINE) + call self % append( ';' // NEWLINE) end subroutine endEntry @@ -227,11 +196,6 @@ subroutine startArray(self, shape) integer(shortInt),dimension(:),intent(in) :: shape character(100), parameter :: Here ='startArray (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state /= IN_ENTRY) then - ! call fatalError(Here,'Cannot finish entry in block or array') - ! end if - ! Update state self % state = IN_ARRAY @@ -240,10 +204,10 @@ subroutine startArray(self, shape) ! Write start of array if( size(self % shapeBuffer) == 1) then - call self % output % append('[ ') + call self % append('[ ') else - call self % output % append('reshape([ ') + call self % append('reshape([ ') end if @@ -259,24 +223,19 @@ subroutine endArray(self) integer(shortInt) :: i character(100), parameter :: Here ='endArray (asciiMATLAB_class.f90)' - ! Check if state support acction - ! if (self % state /= IN_ARRAY) then - ! call fatalError(Here,'Cannot finish array inside an entry') - ! end if - ! Update state self % state = IN_ENTRY ! Write end of array - call self % output % cut(1) - call self % output % append(']') + call self % cut(1) + call self % append(']') ! Finish reshape function for higher rank arrays if( size(self % shapeBuffer) > 1) then do i=1,size(self % shapeBuffer) - call self % output % append(',' // numToChar(self % shapeBuffer(i))) + call self % append(',' // numToChar(self % shapeBuffer(i))) end do - call self % output % append(')') + call self % append(')') end if end subroutine endArray @@ -292,13 +251,11 @@ subroutine printNum(self, val) character(100), parameter :: Here ='printNum (asciiMATLAB_class.f90)' if(self % state == IN_ARRAY) then - call self % output % append(val//',') + call self % append(val//',') else if( self % state == IN_ENTRY) then - call self % output % append(val) + call self % append(val) - ! else - ! call fatalError(Here,'Cannot print number directly into block') end if end subroutine printNum @@ -314,13 +271,11 @@ subroutine printChar(self, val) character(100), parameter :: Here ='printChar (asciiMATLAB_class.f90)' if(self % state == IN_ARRAY) then - call self % output % append(BRAKET_L // APOS // val // APOS // BRAKET_R // ",") + call self % append(BRAKET_L // APOS // val // APOS // BRAKET_R // ",") else if( self % state == IN_ENTRY) then - call self % output % append(BRAKET_L // APOS // val // APOS // BRAKET_R ) + call self % append(BRAKET_L // APOS // val // APOS // BRAKET_R ) - ! else - ! call fatalError(Here,'Cannot print number directly into block') end if end subroutine printChar diff --git a/UserInterface/fileOutput/asciiOutput_inter.f90 b/UserInterface/fileOutput/asciiOutput_inter.f90 index 0b491d921..e60a7b8b2 100644 --- a/UserInterface/fileOutput/asciiOutput_inter.f90 +++ b/UserInterface/fileOutput/asciiOutput_inter.f90 @@ -1,6 +1,7 @@ module asciiOutput_inter use numPrecision + use delayedStream_class, only : delayedStream implicit none private @@ -10,7 +11,7 @@ module asciiOutput_inter !! !! Allows stream-like output to multiple formats (e.g. MATLAB, CSV, JSON) !! - !! They recive the sequence of calls from the `outputFile` and conver them to the output file. + !! They receive the sequence of calls from the `outputFile` and convert them to the output file. !! !! Printers do no error checking. It is responsibility of the `outputFile` to ensure that !! the sequence is correct. @@ -21,7 +22,7 @@ module asciiOutput_inter !! !! Interface: !! init -> Initialise - !! extension -> Return file extension approperiate for the format + !! extension -> Return file extension appropriate for the format !! writeToFile -> Print the output to the provided unit !! startBlock -> Start new block !! endBlock -> End a block @@ -32,10 +33,16 @@ module asciiOutput_inter !! type, public,abstract :: asciiOutput private + type(delayedStream) :: stream contains + procedure :: append + procedure :: cut + procedure :: peek + procedure :: close + procedure :: setUnit procedure(init), deferred :: init + procedure(init), deferred :: endFile procedure(extension), deferred :: extension - procedure(writeToFile),deferred :: writeToFile procedure(startBlock),deferred :: startBlock procedure(endBlock),deferred :: endBlock procedure(startEntry),deferred :: startEntry @@ -61,7 +68,21 @@ subroutine init(self) end subroutine init !! - !! Return approperiate extension for the file + !! End the file + !! + !! Format may require to print some characters at the end, so the printer + !! needs to be notified that no more data will be provided. + !! + !! Args: + !! None + !! + subroutine endFile(self) + import :: asciiOutput + class(asciiOutput), intent(inout) :: self + end subroutine endFile + + !! + !! Return appropriate extension for the file !! !! Must be without any "." Thus "exe" instead of ".exe"! !! @@ -74,19 +95,6 @@ pure function extension(self) result(str) character(:), allocatable :: str end function extension - !! - !! Print the output to the given unit - !! - !! Args: - !! unit [in] -> Unit number of the output - !! - subroutine writeToFile(self, unit) - import :: asciiOutput, & - shortInt - class(asciiOutput), intent(inout) :: self - integer(shortInt), intent(in) :: unit - end subroutine writeToFile - !! !! Change state to writing new block with "name" !! @@ -111,7 +119,7 @@ end subroutine endBlock !! !! Change state to writing a new entry !! - !! Can recive single value or array next + !! Can receive single value or array next !! !! Args: !! name [in] -> name of the entry @@ -134,7 +142,7 @@ end subroutine endEntry !! !! Start writing array with the given shape !! - !! Name should alrady be provided by "startEntry" + !! Name should already be provided by "startEntry" !! !! Args: !! shape [in] -> Shape of the array (in column-major order) @@ -183,4 +191,77 @@ subroutine printChar(self, val) end subroutine printChar end interface +contains + + + !! + !! Write data to the output + !! + !! Args: + !! text [in] -> Text to be written + !! + subroutine append(self, text) + class(asciiOutput), intent(inout) :: self + character(*), intent(in) :: text + + call self % stream % write(text) + + end subroutine append + + !! + !! Remove n characters from the output + !! + !! Args: + !! n [in] -> Number of characters to be removed + !! + subroutine cut(self, n) + class(asciiOutput), intent(inout) :: self + integer, intent(in) :: n + + call self % stream % cut(n) + + end subroutine cut + + !! + !! Get the nth last character from the output + !! + !! Args: + !! n [in] -> How many characters back to peak + !! + function peek(self, n) result(c) + class(asciiOutput), intent(inout) :: self + integer, intent(in) :: n + character(1) :: c + + c = self % stream % peek(n) + + end function peek + + !! + !! Close the stream + !! + !! Signals to the output that no more data will be provided, so the + !! closing characters can be printed. In addition flushes remaining characters + !! in the buffer to the output. + !! + subroutine close(self) + class(asciiOutput), intent(inout) :: self + + call self % endFile() + call self % stream % close() + + end subroutine close + + !! + !! Set the unit to write to + !! + subroutine setUnit(self, unit) + class(asciiOutput), intent(inout) :: self + integer(shortInt), intent(in) :: unit + + call self % stream % setUnit(unit) + + end subroutine setUnit + + end module asciiOutput_inter diff --git a/UserInterface/fileOutput/delayedStream_class.f90 b/UserInterface/fileOutput/delayedStream_class.f90 new file mode 100644 index 000000000..86869fec0 --- /dev/null +++ b/UserInterface/fileOutput/delayedStream_class.f90 @@ -0,0 +1,185 @@ +module delayedStream_class + + use numPrecision + use errors_mod, only : fatalError + + implicit none + private + + integer(shortInt), parameter :: MAX_DELAY = 1024 + integer(shortInt), parameter :: MIN_DELAY = 16 + + !! + !! Writes characters to a file with a delay to allow for backtracking. + !! + !! It is required to make writing of output format in stream-like fashion + !! easier to implement. For example, in JSON, it allows to backtrack and + !! remove trailing comma when closing a JSON dictionary. + !! + !! NOTE: We need to be careful to check if the file is open before writing + !! to it. Otherwise, we will get an error. + !! + !! Private Member: + !! unit -> Unit to which the output file is connected + !! bufferPos -> Position in the buffer + !! buffer -> Buffer of characters + !! + !! Interface: + !! write -> Write some characters to the buffer/file + !! cut -> Remove last N characters from the buffer + !! peek -> Inspect the Nth last character in the buffer + !! close -> Write all characters remaining in the buffer to the file + !! setUnit -> Set the unit to which the output file is connected + !! + type, public :: delayedStream + private + integer(shortInt) :: unit = -8 + integer(shortInt) :: bufferPos = 0 + character(MAX_DELAY) :: buffer + contains + procedure :: write + procedure :: cut + procedure :: peek + procedure :: close + procedure :: setUnit + + procedure, private :: flush + end type delayedStream + +contains + + !! + !! Write to a delayed stream + !! + !! Args: + !! text [in] -> text to write + !! + subroutine write(self, text) + class(delayedStream), intent(inout) :: self + character(*), intent(in) :: text + integer(shortInt) :: i + + ! Write to buffer + do i = 1, len(text) + self % bufferPos = self % bufferPos + 1 + self % buffer(self % bufferPos:self % bufferPos) = text(i:i) + if (self % bufferPos == MAX_DELAY) then + call flush(self) + end if + end do + + end subroutine write + + !! + !! Remove last characters from a delayed stream + !! + !! Args: + !! n [in] -> number of characters to remove. Must be less than the MIN_DELAY. + !! If n is larger than the number of characters in buffer, everything is removed. + !! + subroutine cut(self, n) + class(delayedStream), intent(inout) :: self + integer(shortInt), intent(in) :: n + character(100), parameter :: HERE = "cut (delayedStream_class.f90)" + + if (n > MIN_DELAY) then + call fatalError(HERE, "n must be less than MIN_DELAY") + end if + + self % bufferPos = max(0, self % bufferPos - n) + + end subroutine cut + + !! + !! Inspect the nth last character in the buffer + !! + !! Args: + !! n [in] -> The number of characters to peek back. + !! + function peek(self, n) result(c) + class(delayedStream), intent(inout) :: self + integer(shortInt), intent(in) :: n + character(n) :: c + character(100), parameter :: HERE = "peek (delayedStream_class.f90)" + + if (n > self % bufferPos) then + call fatalError(HERE, "n is larger than the buffer") + else + c = self % buffer(self % bufferPos - n + 1: self % bufferPos - n + 1) + end if + + end function peek + + !! + !! Flush contents of the buffer + !! + subroutine flush(self) + class(delayedStream), intent(inout) :: self + integer(shortInt) :: to_file + logical(defBool) :: isOpen + + ! Do not write anything if the file is not opened + inquire(unit=self % unit, opened=isOpen) + if (.not. isOpen) then + return + end if + + ! Do nothing if the buffer is below minimum length + if (self % bufferPos <= MIN_DELAY) then + return + end if + to_file = self % bufferPos - MIN_DELAY + + ! Write to file only if it is open + if (isOpen) then + write(self % unit, "(A)", advance="no") self % buffer(1 : to_file) + end if + + self % buffer(1:MIN_DELAY) = self % buffer(to_file + 1 : self % bufferPos) + + self % bufferPos = MIN_DELAY + + end subroutine flush + + !! + !! Close the delayed stream + !! + subroutine close(self) + class(delayedStream), intent(inout) :: self + logical(defBool) :: isOpen + + ! Do not write anything if the file is not opened + inquire(unit=self % unit, opened=isOpen) + if (.not. isOpen) then + return + end if + + ! Does nothing if empty + if (self % bufferPos == 0) then + return + end if + + ! Write the remaining buffer + ! Write to file only if it is open + if (isOpen) then + write(self % unit, "(A)", advance="no") self % buffer(1:self % bufferPos) + end if + self % bufferPos = 0 + + end subroutine close + + !! + !! Set the unit to which the output file is connected + !! + !! Args: + !! unit [in] -> Unit to which the output file is connected. May be opened or not. + !! Both are fine. If it is not opened, nothing will be written to it. + subroutine setUnit(self, unit) + class(delayedStream), intent(inout) :: self + integer(shortInt), intent(in) :: unit + + self % unit = unit + + end subroutine setUnit + +end module delayedStream_class diff --git a/UserInterface/fileOutput/dummyPrinter_class.f90 b/UserInterface/fileOutput/dummyPrinter_class.f90 index 890883d4d..de2f84442 100644 --- a/UserInterface/fileOutput/dummyPrinter_class.f90 +++ b/UserInterface/fileOutput/dummyPrinter_class.f90 @@ -23,8 +23,8 @@ module dummyPrinter_class contains procedure :: init + procedure :: endFile procedure :: extension - procedure :: writeToFile procedure :: startBlock procedure :: endBlock @@ -52,29 +52,29 @@ subroutine init(self) end subroutine init !! - !! Return approperiate extension for the file + !! End the file !! !! See asciiOutput_inter for details !! - pure function extension(self) result(str) - class(dummyPrinter), intent(in) :: self - character(:), allocatable :: str + subroutine endFile(self) + class(dummyPrinter), intent(inout) :: self - str = '' + ! Nothing to do - end function extension + end subroutine endFile !! - !! Print the output to the given unit + !! Return approperiate extension for the file !! !! See asciiOutput_inter for details !! - subroutine writeToFile(self, unit) - class(dummyPrinter), intent(inout) :: self - integer(shortInt), intent(in) :: unit + pure function extension(self) result(str) + class(dummyPrinter), intent(in) :: self + character(:), allocatable :: str + str = '' - end subroutine writeToFile + end function extension !! !! Change state to writing new block with "name" diff --git a/UserInterface/fileOutput/outputFile_class.f90 b/UserInterface/fileOutput/outputFile_class.f90 index 6585580e0..48fb2011d 100644 --- a/UserInterface/fileOutput/outputFile_class.f90 +++ b/UserInterface/fileOutput/outputFile_class.f90 @@ -1,7 +1,8 @@ module outputFile_class use numPrecision - use genericProcedures, only : fatalError, numToChar + use errors_mod, only : fatalError + use genericProcedures, only : numToChar use stack_class, only : stackChar use charTape_class, only : charTape use charMap_class, only : charMap @@ -28,10 +29,10 @@ module outputFile_class !! !! Stack to store the dictionaries with occupied entries !! - !! Note that blocks in outputFile are enumarated from 0 + !! Note that blocks in outputFile are enumerated from 0 !! So index in stack is idx = blockLevel + 1 !! - !! Imlementation must be rebust to avoid segmentation errors if outut is used + !! Implementation must be robust to avoid segmentation errors if output is used !! in incorrect sequence (e.g. closing more blocks then were open) !! !! Interface @@ -69,7 +70,7 @@ module outputFile_class !! is called root !! -> Each entry in a block (including sub-blocks) must have a unique name !! -> Each entry is either a raw "Name - Value" pair or an array - !! -> Arrays can have unlimited rank (dimension), but values of the enteries must be given in a + !! -> Arrays can have unlimited rank (dimension), but values of the entries must be given in a !! sequence in a COLUMN-MAJOR ORDER !! -> Arrays can store RESULTS, REALS, INTS or CHARACTERS. !! -> RESULT is a pair of reals. MEAN and STANDARD DEVIATION @@ -81,14 +82,14 @@ module outputFile_class !! !! Handling errors: !! If `fatalError` is true (default) then on wrong sequence of calls fatalError will be raised. - !! Otherwise the flag `noErrors` will be set to false and error massege saved in the `errorLog` + !! Otherwise the flag `noErrors` will be set to false and error massage saved in the `errorLog` !! !! Repeated Names: - !! If name in a block is repeted bahviour is governed by `logNameReuse` function. If `fatalError` - !! is true, the nwarning message is produced. Otherwise error is logged in `errorLog` and + !! If name in a block is repeated behaviour is governed by `logNameReuse` function. If `fatalError` + !! is true, the warning message is produced. Otherwise error is logged in `errorLog` and !! `noErrors` flag set to FALSE. !! - !! Private Memebers: + !! Private Members: !! output -> Output printer that determines the format !! type -> Type of the output printer !! fatalErrors -> True (default) makes all errors fatal @@ -107,8 +108,7 @@ module outputFile_class !! !! Interface: !! init -> Initialise the output file with format - !! writeToFile -> Write the output to a file - !! writeToConsole -> Write the output to the console + !! finalise -> Finalise the output file. Ensures all data is written. !! reset -> Discard already written output. !! isValid -> Return true if no errors were raised since initialisation !! getErrorLog -> Get the error log @@ -124,8 +124,10 @@ module outputFile_class !! type, public :: outputFile private - class(asciiOutput),allocatable :: output - character(nameLen) :: type + class(asciiOutput), allocatable :: output + character(nameLen) :: type + character(:), allocatable :: outputFileName + integer(shortInt) :: outputUnit = -97875674 ! Hopefully this does not exist for any internal units ! Error handling settings logical(defBool) :: fatalErrors = .true. ! Enable fatalErrors on wrong logic @@ -141,7 +143,7 @@ module outputFile_class character(nameLen) :: current_array_name = '' type(stackChar) :: block_name_stack - ! Buffors + ! Buffers integer(shortInt),dimension(:), allocatable :: shapeBuffer type(charMapStack) :: usedNames @@ -151,9 +153,8 @@ module outputFile_class contains procedure :: init - procedure :: writeToFile - procedure :: writeToConsole procedure :: reset + final :: finalisation ! Error Handling procedure, private :: logError @@ -210,17 +211,56 @@ module outputFile_class !! !! Args: !! type [in] -> Character to specify the format of the output file - !! fatalErrors [in] -> Optional. Set on/off fatalErrors oninvalid use of output file. + !! fatalErrors [in] -> Optional. Set on/off fatalErrors on invalid use of output file. !! Default is ON. - !! - subroutine init(self, type, fatalErrors) - class(outputFile), intent(inout) :: self - character(*),intent(in) :: type - logical(defBool),optional,intent(in) :: fatalErrors + !! filename [in] -> Optional. Set the name of the output file. If none is given no output is + !! written. Must be given without extension. It will be added automatically + !! based on the format of the output file. + !! + subroutine init(self, type, fatalErrors, filename) + class(outputFile), intent(inout) :: self + character(*), intent(in) :: type + logical(defBool), optional,intent(in) :: fatalErrors + character(*), optional, intent(in) :: filename + integer(shortInt) :: error + character(99) :: errorMsg + logical(defBool) :: isOpen + character(100), parameter :: Here = "init (outputFile_class.f90)" self % type = type allocate( self % output, source = new_asciiOutput(self % type)) + ! Open file if given + ! Otherwise make sure that the unit is not opened + if (present(filename)) then + ! Because GFotran hates implicit allocation + ! It causes buggy warnings to appear in compilation + allocate(self % outputFileName, source = trim(filename) // "." // self % output % extension()) + + open(newunit = self % outputUnit, & + file = self % outputFileName, & + status = "replace", & + action = "write", & + access = "stream", & + form = "formatted",& + iostat = error, & + iomsg = errorMsg) + + ! Catch error + if (error /= 0 ) call fatalError(Here, errorMsg) + else + ! Set unit to some unused, unopened unit + ! delayedStream will not write anything if given an unopened file + do + inquire(self % outputUnit, opened = isOpen) + if (.not. isOpen) exit + self % outputUnit = self % outputUnit + 1 + end do + end if + + ! Set output unit + call self % output % setUnit(self % outputUnit) + ! Initialise name stack call self % usedNames % init() call self % usedNames % push() @@ -233,58 +273,36 @@ subroutine init(self, type, fatalErrors) end subroutine init !! - !! Dump collected results to file - !! - !! Writes the output stored in memory to a file - !! - !! Args: - !! file [in] -> Path to the output file. Must be given WITHOUT extension - !! Approperiate extension will recived from the printer + !! Upon the destruction of the object, write the output to file !! - !! Errors: - !! fatalError if there is a problem opening the file. - !! - subroutine writeToFile(self, file) - class(outPutFile), intent(inout) :: self - character(*), intent(in) :: file - integer(shortInt) :: unitNum - integer(shortInt) :: error - character(:), allocatable :: file_loc - character(99) :: errorMsg - character(100), parameter :: Here = 'writeToFile (outputFile_class.f90)' - - file_loc = trim(file) // "." // self % output % extension() - - ! Open file to write - open ( newunit = unitNum, & - file = file_loc, & - action = 'write', & - iostat = error , & - iomsg = errorMsg ) - - ! Catch error - if (error /= 0 ) call fatalError(Here, errorMsg) - - ! Write to file - call self % output % writeToFile(unitNum) - - ! Close file - close(unitNum) + !! NOTE: + !! Being marked `final`(a deconstructor) this subroutine will be normally implicitly called by + !! by the compiler whenever a variable of `type(outputFile)` goes out of scope (e.g. when program + !! execution reaches `end subroutine`, `end function` or `end block` statement. However, the rule in + !! Fortran is that the decostructors are NOT called at the end of the program `end program`. In that case + !! if one uses the output file as a module variable or defines it inside the program block, the `finalisation` + !! must be called explicitly to ensure or data from the buffer is written to the disk + !! + !! + subroutine finalisation(self) + type(outputFile), intent(inout) :: self + logical(defBool) :: isOpen + + !! Deallocate the printer + !! Needs to be done before closing the file + !! So remaining contents of the buffer are flushed + if (allocated(self % output)) then + call self % output % close() + deallocate(self % output) + end if - end subroutine writeToFile + inquire(unit = self % outputUnit, opened = isOpen) - !! - !! Print output file to the console (default output) - !! - !! Args: - !! None - !! - subroutine writeToConsole(self) - class(outputFile), intent(inout) :: self - - call self % output % writeToFile(OUTPUT_UNIT) + if (isOpen) then + close(self % outputUnit) + end if - end subroutine writeToConsole + end subroutine finalisation !! !! Discards all results and sets output file to an initialised state @@ -296,6 +314,8 @@ subroutine reset(self) deallocate(self % output) allocate( self % output, source = new_asciiOutput(self % type)) + if (allocated(self % outputFileName)) deallocate(self % outputFileName) + ! Clean error log an reset flag call self % errorLog % clean() self % noErrors = .true. @@ -341,10 +361,10 @@ subroutine logError(self, where, what) end subroutine logError !! - !! Deal with an event when entry name is resued in a block + !! Deal with an event when entry name is reused in a block !! !! Currently is fatalError is TRUE -> Prints warning to the screen - !! If fatalError is False markes error flagto true and appends a log + !! If fatalError is False marks error flag to true and appends the log !! !! Args: !! name [in] -> Name that has been reused @@ -403,7 +423,7 @@ end function getErrorLog !! !! Errors: !! Calls logError if: - !! -Called when writting an array + !! -Called when writing an array !! Calls logNameReuse if name is not unique !! subroutine startBlock(self, name) @@ -422,7 +442,7 @@ subroutine startBlock(self, name) if ( self % arrayType /= NOT_ARRAY) then call self % logError(Here,'In block: '// self % current_block_name // & ' Cannot begin new block: '//name// & - ' Becouse is writing array: '// self % current_array_name) + ' Because is writing array: '// self % current_array_name) end if ! Update state @@ -443,7 +463,7 @@ end subroutine startBlock !! !! Errors: !! Calls logError if: - !! -Called when writting an array + !! -Called when writing an array !! -Tries to end the root block. !! subroutine endBlock(self) @@ -454,7 +474,7 @@ subroutine endBlock(self) if ( self % arrayType /= NOT_ARRAY) then call self % logError(Here,'In block: '// self % current_block_name // & ' Cannot exit from block: '// & - ' Becouse is writing array: '// self % current_array_name) + ' Because is writing array: '// self % current_array_name) end if ! Check that is not in root block @@ -485,7 +505,7 @@ end subroutine endBlock !! !! Errors: !! Calls logError if: - !! -Called when writting another array + !! -Called when writing another array !! -Shape is invalid !! Calls logNameReuse if name is not unique !! @@ -505,14 +525,14 @@ subroutine startArray(self, name, shape) if ( self % arrayType /= NOT_ARRAY) then call self % logError(Here,'In block: '// self % current_block_name // & ' Cannot start new array: '//name// & - ' Becouse is writing array: '// self % current_array_name) + ' Because is writing array: '// self % current_array_name) end if ! Check whether shape is degenerate if (size(shape) == 0 .or. product(shape) == 0) then call self % logError(Here,'In block: '// self % current_block_name // & ' Cannot start new array: '//name// & - ' Becouse the shape is degenerate') + ' Because the shape is degenerate') end if ! Load shape information into buffer @@ -524,8 +544,8 @@ subroutine startArray(self, name, shape) ! Update State -> Array is present but has undefined type self % arrayType = UNDEF_ARRAY - ! Print begining of new entry - ! We don't know the shape of array becouse if has results its rank will increase + ! Print beginning of new entry + ! We don't know the shape of array because if has results its rank will increase call self % output % startEntry(name) end subroutine startArray @@ -548,7 +568,7 @@ subroutine endArray(self) if(self % arrayTop /= self % arrayLimit) then call self % logError(Here, 'Cannot close array: ' // trim(self % current_array_name) // & ' in block: ' // trim(self % current_block_name) // & - ' Becouse only: '// numToChar(self % arrayTop) // & + ' Because only: '// numToChar(self % arrayTop) // & ' of '// numToChar(self % arrayLimit)//' were provided') end if @@ -615,10 +635,10 @@ end subroutine addResult_scalar !! Add results to the array !! !! Args: - !! val [in] -> Array of mean valuse + !! val [in] -> Array of mean values !! std [in] -> Associated standard deviations !! - !! Erorrs: + !! Errors: !! Calls logError if: !! -Called before an array is opened !! -Currently opened array has different type @@ -632,7 +652,7 @@ subroutine addResult_rank1(self, val, std) character(100),parameter :: Here ='addResult_rank1 (outputFile_class.f90)' N = size(val) - if (N /= size(std)) call self % logError(Here,'val and std have diffrent size.') + if (N /= size(std)) call self % logError(Here, 'val and std have different size.') ! Add all individual entries do i = 1, N @@ -1203,7 +1223,7 @@ end subroutine push_charMapStack !! None !! !! Errors: - !! Poping empty stack does nothing! + !! Popping empty stack does nothing! !! subroutine pop_charMapStack(self) class(charMapStack), intent(inout) :: self @@ -1218,7 +1238,7 @@ subroutine pop_charMapStack(self) end subroutine pop_charMapStack !! - !! Add new entery at the given block level + !! Add new entry at the given block level !! !! Args: !! name [in] -> Entry name @@ -1246,7 +1266,7 @@ end subroutine add_charMapStack !! b_lvl [in] -> Block level (idx = b_lvl + 1) !! !! Result: - !! True if it is. False otherwsie + !! True if it is. False otherwise !! !! Errors: !! Returns FALSE is b_lvl is out of range diff --git a/Visualisation/visualiser_class.f90 b/Visualisation/visualiser_class.f90 index 9641f7758..a0cdcd768 100644 --- a/Visualisation/visualiser_class.f90 +++ b/Visualisation/visualiser_class.f90 @@ -1,394 +1,398 @@ -module visualiser_class - - use numPrecision - use universalVariables - use genericProcedures, only : fatalError, numToChar - use hashFunctions_func, only : knuthHash - use imgBmp_func, only : imgBmp_toFile - use commandLineUI, only : getInputFile - use dictionary_class, only : dictionary - use geometry_inter, only : geometry - use outputVTK_class - - implicit none - private - - !! - !! Object responsible for controlling visualisation - !! - !! Object that creates images relating to SCONE geometries - !! Should be extensible for adding different visualisation methods - !! Recieves and generates data for visualisation - !! Requires a dictionary input which specifies the procedures to call - !! Presently supports: VTK voxel mesh creation - !! - !! Private members: - !! name -> name to be used for generating output files (corresponding to input) - !! geom -> pointer to geometry - !! vizDict -> dictionary containing visualisations to be generated - !! - !! Interface: - !! init -> initialises visualiser - !! makeViz -> constructs requested visualisations - !! kill -> cleans up visualiser - !! - !! Sample dictionary input: - !! viz{ - !! vizDict1{ } - !! #vizDict2{ }# - !! } - !! - !! NOTE: For details regarding contents of the vizDict dictionaries see the documentation - !! of 'makeVTK' and 'makeBmpImg' functions - !! - type, public :: visualiser - character(nameLen), private :: name - class(geometry), pointer, private :: geom => null() - type(dictionary), private :: vizDict - contains - procedure :: init - procedure :: makeViz - procedure :: kill - procedure, private :: makeVTK - procedure, private :: makeBmpImg - end type - -contains - - !! - !! Initialises visualiser - !! - !! Provides visualiser with filename for output, - !! geometry information, and the dictionary decribing - !! what is to be plotted - !! - !! Args: - !! geom [inout] -> pointer to the geometry - !! vizDict[in] -> dictionary containing what is to be visualised - !! - !! Result: - !! Initialised visualiser - !! - subroutine init(self, geom, vizDict) - class(visualiser), intent(inout) :: self - class(geometry), pointer, intent(inout) :: geom - class(dictionary), intent(in) :: vizDict - character(:), allocatable :: string - - ! Obtain file name - call getInputFile(string) - self % name = string - - ! Point to geometry - self % geom => geom - - ! Store visualisation dictionary - self % vizDict = vizDict - - end subroutine init - - !! - !! Generate all visualisations specified by vizDict - !! - !! Proceed through all dictionaries contained within vizDict - !! and perform all corresponding visualisations - !! - !! Result: - !! Visualisation outputs corresponding to dictionary contents - !! - !! Errors: - !! Returns an error if an unrecognised visualisation is requested - !! - subroutine makeViz(self) - class(visualiser), intent(inout) :: self - class(dictionary), pointer :: tempDict - character(nameLen),dimension(:), allocatable :: keysArr - integer(shortInt) :: i - character(nameLen) :: type - character(nameLen) :: here ='makeViz (visualiser_class.f90)' - - ! Loop through each sub-dictionary and generate visualisation - ! (if the visualisation method is available) - call self % vizDict % keys(keysArr,'dict') - - do i=1,size(keysArr) - tempDict => self % vizDict % getDictPtr(keysArr(i)) - call tempDict % get(type,'type') - select case(type) - case('vtk') - call self % makeVTK(tempDict) - - case('bmp') - call self % makeBmpImg(tempDict) - - case default - call fatalError(here, 'Unrecognised visualisation - presently only accept vtk') - - end select - - end do - - end subroutine makeViz - - !! - !! Generate a VTK output - !! - !! Creates the VTK file corresponding to the contents of dict - !! - !! Args: - !! dict [in] -> dictionary containing description of VTK file to be made - !! - !! Sample input dictionary: - !! VTK { - !! type vtk; - !! corner (-1.0 -1.0 -1.0); // lower corner of the plot volume - !! width (2.0 2.0 2.0); // width in each direction - !! vox (300 300 300); // Resolution in each direction - !! #what uniqueId;# // Plot target. 'material' or 'uniqueId'. Default: 'material' - !! } - !! - !! TODO: VTK output is placed in a input filename appended by '.vtk' extension. - !! This prevents multiple VTK visualistions (due to overriding). Might also become - !! weird for input files with extension e.g. 'input.dat'. - !! DEMAND USER TO GIVE OUTPUT NAME - !! - subroutine makeVTK(self, dict) - class(visualiser), intent(inout) :: self - class(dictionary), intent(in) :: dict - type(outputVTK) :: vtk - integer(shortInt), dimension(:,:,:), allocatable:: voxelMat - real(defReal), dimension(:), allocatable :: corner ! corner of the mesh - real(defReal), dimension(:), allocatable :: center ! center of the mesh - real(defReal), dimension(:), allocatable :: width ! corner of the mesh - integer(shortInt), dimension(:), allocatable :: nVox ! number of mesh voxels - character(nameLen) :: what - character(nameLen) :: here ='makeVTK (visualiser_class.f90)' - - call vtk % init(dict) - - ! Identify whether plotting 'material' or 'cellID' - call dict % getOrDefault(what, 'what', 'material') - - ! Obtain geometry data - call dict % get(corner, 'corner') - call dict % get(width, 'width') - center = corner + width/TWO - call dict % get(nVox, 'vox') - - if (size(corner) /= 3) then - call fatalError(here,'Voxel plot requires corner to have 3 values') - endif - if (size(width) /= 3) then - call fatalError(here,'Voxel plot requires width to have 3 values') - endif - if (size(nVox) /= 3) then - call fatalError(here,'Voxel plot requires vox to have 3 values') - endif - allocate(voxelMat(nVox(1), nVox(2), nVox(3))) - - ! Have geometry obtain data - call self % geom % voxelPlot(voxelMat, center, what, width) - - ! In principle, can add multiple data sets to VTK - not done here yet - ! VTK data set will use 'what' variable as a name - call vtk % addData(voxelMat, what) - call vtk % output(self % name) - call vtk % kill() - - end subroutine makeVTK - - !! - !! Generate a BMP slice image of the geometry - !! - !! Args: - !! dict [in] -> Dictionary with settings - !! - !! Sample dictionary input: - !! bmp_img { - !! type bmp; - !! #what uniqueID;# // Target of the plot. 'uniqueId' or 'material'. Default: 'material' - !! output img; // Name of output file without extension - !! centre (0.0 0.0 0.0); // Coordinates of the centre of the plot - !! axis x; // Must be 'x', 'y' or 'z' - !! res (300 300); // Resolution of the image - !! #width (1.0 2.0);# // Width of the plot from the centre - !! } - !! - !! NOTE: If 'width' is not given, the plot will extend to the bounds of the geometry. - !! This may result in the provided centre beeing moved to the center of the geoemtry in the - !! plot plane. However, the position on the plot axis will be unchanged. - !! - subroutine makeBmpImg(self, dict) - class(visualiser), intent(inout) :: self - class(dictionary), intent(in) :: dict - real(defReal), dimension(3) :: centre - real(defReal), dimension(2) :: width - character(1) :: dir - character(nameLen) :: tempChar - logical(defBool) :: useWidth - character(nameLen) :: what, outputFile - real(defReal), dimension(:), allocatable :: temp - integer(shortInt), dimension(:), allocatable :: tempInt - integer(shortInt), dimension(:,:), allocatable :: img - character(100), parameter :: Here = 'makeBmpImg (visualiser_class.f90)' - - ! Get plot parameters - - ! Identify whether plotting 'material' or 'cellID' - call dict % getOrDefault(what, 'what', 'material') - - ! Get name of the output file - call dict % get(outputFile, 'output') - outputFile = trim(outputFile) // '.bmp' - - ! Central point - call dict % get(temp, 'centre') - - if (size(temp) /= 3) then - call fatalError(Here, "'center' must have size 3. Has: "//numToChar(size(temp))) - end if - - centre = temp - - ! Axis - call dict % get(tempChar, 'axis') - - if (len_trim(tempChar) /= 1) then - call fatalError(Here, "'axis' must be x,y or z. Not: "//tempChar) - end if - - dir = tempChar(1:1) - - ! Resolution - call dict % get(tempInt, 'res') - - if (size(tempInt) /= 2) then - call fatalError(Here, "'res' must have size 2. Has: "//numToChar(size(tempInt))) - else if (any(tempInt <= 0)) then - call fatalError(Here, "Resolution must be +ve. There is 0 or -ve entry!") - end if - - allocate(img(tempInt(1), tempInt(2))) - - ! Optional width - useWidth = dict % isPresent('width') - if (useWidth) then - call dict % get(temp, 'width') - - ! Check for errors - if (size(temp) /= 2) then - call fatalError(Here, "'width' must have size 2. Has: "//numToChar((size(temp)))) - else if (any(temp <= ZERO)) then - call fatalError(Here, "'width' must be +ve. It isn't.") - end if - - width = temp - - end if - - ! Get plot - if (useWidth) then - call self % geom % slicePlot(img, centre, dir, what, width) - else - call self % geom % slicePlot(img, centre, dir, what) - end if - - ! Translate to an image - select case (what) - case ('material') - img = materialColor(img) - - case ('uniqueID') - img = uniqueIDColor(img) - - case default - call fatalError(Here, "Invalid request for plot target. Must be 'material' or 'uniqueID'& - & is: "//what) - end select - - ! Print image - call imgBmp_toFile(img, outputFile) - - end subroutine makeBmpImg - - !! - !! Terminates visualiser - !! - !! Cleans up remnants of visualiser once it is no longer needed - !! - !! Result: - !! An empty visualiser object - !! - subroutine kill(self) - class(visualiser), intent(inout) :: self - - self % name ='' - self % geom => null() - call self % vizDict % kill() - - end subroutine kill - - - !! - !! Convert matIdx to a 24bit color - !! - !! Special materials are associeted with special colors: - !! OUTSIDE_MAT -> white (#ffffff) - !! VOID_MAT -> black (#000000) - !! UNDEF_MAT -> green (#00ff00) - !! - !! Args: - !! matIdx [in] -> Value of the material index - !! - !! Result: - !! A 24-bit color specifing the material - !! - elemental function materialColor(matIdx) result(color) - integer(shortInt), intent(in) :: matIdx - integer(shortInt) :: color - integer(shortInt), parameter :: COL_OUTSIDE = int(z'ffffff', shortInt) - integer(shortInt), parameter :: COL_VOID = int(z'000000', shortInt) - integer(shortInt), parameter :: COL_UNDEF = int(z'00ff00', shortInt) - - select case (matIdx) - case (OUTSIDE_MAT) - color = COL_OUTSIDE - - case (VOID_MAT) - color = COL_VOID - - case (UNDEF_MAT) - color = COL_UNDEF - - case default - color = knuthHash(matIdx, 24) - - end select - - end function materialColor - - !! - !! Convert uniqueID to 24bit color - !! - !! An elemental wrapper over Knuth Hash - !! - !! We use a hash function to scatter colors accross all available. - !! Knuth multiplicative hash is very good at scattering integer - !! sequences e.g. {1, 2, 3...}. Thus, it is ideal for a colormap. - !! - !! Args: - !! uniqueID [in] -> Value of the uniqueID - !! - !! Result: - !! A 24-bit color specifing the uniqueID - !! - elemental function uniqueIDColor(uniqueID) result(color) - integer(shortInt), intent(in) :: uniqueID - integer(shortInt) :: color - - color = knuthHash(uniqueID, 24) - - end function uniqueIDColor - - -end module visualiser_class +module visualiser_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, numToChar + use hashFunctions_func, only : knuthHash, FNV_1 + use imgBmp_func, only : imgBmp_toFile + use commandLineUI, only : getInputFile + use dictionary_class, only : dictionary + use geometry_inter, only : geometry + use materialMenu_mod, only : mm_colourMap => colourMap + use outputVTK_class + + implicit none + private + + !! + !! Object responsible for controlling visualisation + !! + !! Object that creates images relating to SCONE geometries + !! Should be extensible for adding different visualisation methods + !! Receives and generates data for visualisation + !! Requires a dictionary input which specifies the procedures to call + !! Presently supports: VTK voxel mesh creation + !! + !! Private members: + !! name -> name to be used for generating output files (corresponding to input) + !! geom -> pointer to geometry + !! vizDict -> dictionary containing visualisations to be generated + !! + !! Interface: + !! init -> initialises visualiser + !! makeViz -> constructs requested visualisations + !! kill -> cleans up visualiser + !! + !! Sample dictionary input: + !! viz{ + !! vizDict1{ } + !! #vizDict2{ }# + !! } + !! + !! NOTE: For details regarding contents of the vizDict dictionaries see the documentation + !! of 'makeVTK' and 'makeBmpImg' functions + !! + type, public :: visualiser + character(nameLen), private :: name + class(geometry), pointer, private :: geom => null() + type(dictionary), private :: vizDict + contains + procedure :: init + procedure :: makeViz + procedure :: kill + procedure, private :: makeVTK + procedure, private :: makeBmpImg + end type + +contains + + !! + !! Initialises visualiser + !! + !! Provides visualiser with filename for output, + !! geometry information, and the dictionary describing + !! what is to be plotted + !! + !! Args: + !! geom [inout] -> pointer to the geometry + !! vizDict[in] -> dictionary containing what is to be visualised + !! + !! Result: + !! Initialised visualiser + !! + subroutine init(self, geom, vizDict) + class(visualiser), intent(inout) :: self + class(geometry), pointer, intent(inout) :: geom + class(dictionary), intent(in) :: vizDict + character(:), allocatable :: string + + ! Obtain file name + call getInputFile(string) + self % name = string + + ! Point to geometry + self % geom => geom + + ! Store visualisation dictionary + self % vizDict = vizDict + + end subroutine init + + !! + !! Generate all visualisations specified by vizDict + !! + !! Proceed through all dictionaries contained within vizDict + !! and perform all corresponding visualisations + !! + !! Result: + !! Visualisation outputs corresponding to dictionary contents + !! + !! Errors: + !! Returns an error if an unrecognised visualisation is requested + !! + subroutine makeViz(self) + class(visualiser), intent(inout) :: self + class(dictionary), pointer :: tempDict + character(nameLen),dimension(:), allocatable :: keysArr + integer(shortInt) :: i + character(nameLen) :: type + character(nameLen) :: here ='makeViz (visualiser_class.f90)' + + ! Loop through each sub-dictionary and generate visualisation + ! (if the visualisation method is available) + call self % vizDict % keys(keysArr,'dict') + + do i=1,size(keysArr) + tempDict => self % vizDict % getDictPtr(keysArr(i)) + call tempDict % get(type,'type') + select case(type) + case('vtk') + call self % makeVTK(tempDict) + + case('bmp') + call self % makeBmpImg(tempDict) + + case default + call fatalError(here, 'Unrecognised visualisation - presently only accept vtk') + + end select + + end do + + end subroutine makeViz + + !! + !! Generate a VTK output + !! + !! Creates the VTK file corresponding to the contents of dict + !! + !! Args: + !! dict [in] -> dictionary containing description of VTK file to be made + !! + !! Sample input dictionary: + !! VTK { + !! type vtk; + !! corner (-1.0 -1.0 -1.0); // lower corner of the plot volume + !! width (2.0 2.0 2.0); // width in each direction + !! vox (300 300 300); // Resolution in each direction + !! #what uniqueId;# // Plot target. 'material' or 'uniqueId'. Default: 'material' + !! } + !! + !! TODO: VTK output is placed in a input filename appended by '.vtk' extension. + !! This prevents multiple VTK visualisations (due to overriding). Might also become + !! weird for input files with extension e.g. 'input.dat'. + !! DEMAND USER TO GIVE OUTPUT NAME + !! + subroutine makeVTK(self, dict) + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + type(outputVTK) :: vtk + integer(shortInt), dimension(:,:,:), allocatable:: voxelMat + real(defReal), dimension(:), allocatable :: corner ! corner of the mesh + real(defReal), dimension(:), allocatable :: center ! center of the mesh + real(defReal), dimension(:), allocatable :: width ! corner of the mesh + integer(shortInt), dimension(:), allocatable :: nVox ! number of mesh voxels + character(nameLen) :: what + character(nameLen) :: here ='makeVTK (visualiser_class.f90)' + + call vtk % init(dict) + + ! Identify whether plotting 'material' or 'cellID' + call dict % getOrDefault(what, 'what', 'material') + + ! Obtain geometry data + call dict % get(corner, 'corner') + call dict % get(width, 'width') + center = corner + width/TWO + call dict % get(nVox, 'vox') + + if (size(corner) /= 3) then + call fatalError(here,'Voxel plot requires corner to have 3 values') + endif + if (size(width) /= 3) then + call fatalError(here,'Voxel plot requires width to have 3 values') + endif + if (size(nVox) /= 3) then + call fatalError(here,'Voxel plot requires vox to have 3 values') + endif + allocate(voxelMat(nVox(1), nVox(2), nVox(3))) + + ! Have geometry obtain data + call self % geom % voxelPlot(voxelMat, center, what, width) + + ! In principle, can add multiple data sets to VTK - not done here yet + ! VTK data set will use 'what' variable as a name + call vtk % addData(voxelMat, what) + call vtk % output(self % name) + call vtk % kill() + + end subroutine makeVTK + + !! + !! Generate a BMP slice image of the geometry + !! + !! Args: + !! dict [in] -> Dictionary with settings + !! + !! Sample dictionary input: + !! bmp_img { + !! type bmp; + !! #what uniqueID;# // Target of the plot. 'uniqueId' or 'material'. Default: 'material' + !! output img; // Name of output file without extension + !! centre (0.0 0.0 0.0); // Coordinates of the centre of the plot + !! axis x; // Must be 'x', 'y' or 'z' + !! res (300 300); // Resolution of the image + !! #offset 978; # // Parameter to 'randomize' the colour map + !! #width (1.0 2.0);# // Width of the plot from the centre + !! } + !! + !! NOTE: If 'width' is not given, the plot will extend to the bounds of the geometry. + !! This may result in the provided centre being moved to the center of the geometry in the + !! plot plane. However, the position on the plot axis will be unchanged. + !! + subroutine makeBmpImg(self, dict) + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal), dimension(3) :: centre + real(defReal), dimension(2) :: width + character(1) :: dir + character(nameLen) :: tempChar + logical(defBool) :: useWidth + character(nameLen) :: what, outputFile + real(defReal), dimension(:), allocatable :: temp + integer(shortInt), dimension(:), allocatable :: tempInt + integer(shortInt), dimension(:,:), allocatable :: img + integer(shortInt) :: offset + character(10) :: time + character(8) :: date + character(100), parameter :: Here = 'makeBmpImg (visualiser_class.f90)' + + ! Get plot parameters + + ! Identify whether plotting 'material' or 'cellID' + call dict % getOrDefault(what, 'what', 'material') + + ! Get name of the output file + call dict % get(outputFile, 'output') + outputFile = trim(outputFile) // '.bmp' + + ! Central point + call dict % get(temp, 'centre') + + if (size(temp) /= 3) then + call fatalError(Here, "'centre' must have size 3. Has: "//numToChar(size(temp))) + end if + + centre = temp + + ! Axis + call dict % get(tempChar, 'axis') + + if (len_trim(tempChar) /= 1) then + call fatalError(Here, "'axis' must be x,y or z. Not: "//tempChar) + end if + + dir = tempChar(1:1) + + ! Resolution + call dict % get(tempInt, 'res') + + if (size(tempInt) /= 2) then + call fatalError(Here, "'res' must have size 2. Has: "//numToChar(size(tempInt))) + else if (any(tempInt <= 0)) then + call fatalError(Here, "Resolution must be +ve. There is 0 or -ve entry!") + end if + + allocate(img(tempInt(1), tempInt(2))) + + ! Optional width + useWidth = dict % isPresent('width') + if (useWidth) then + call dict % get(temp, 'width') + + ! Check for errors + if (size(temp) /= 2) then + call fatalError(Here, "'width' must have size 2. Has: "//numToChar((size(temp)))) + else if (any(temp <= ZERO)) then + call fatalError(Here, "'width' must be +ve. It isn't.") + end if + + width = temp + + end if + + ! Colourmap offset + ! If not given select randomly + if (dict % isPresent('offset')) then + call dict % get(offset, 'offset') + + else + call date_and_time(date, time) + call FNV_1(date // time, offset) + + end if + + ! Get plot + if (useWidth) then + call self % geom % slicePlot(img, centre, dir, what, width) + else + call self % geom % slicePlot(img, centre, dir, what) + end if + + ! Translate to an image + select case (what) + case ('material') + img = materialColour(img, offset) + + case ('uniqueID') + img = uniqueIDColour(img) + + case default + call fatalError(Here, "Invalid request for plot target. Must be 'material' or 'uniqueID'& + & is: "//what) + end select + + ! Print image + call imgBmp_toFile(img, outputFile) + + end subroutine makeBmpImg + + !! + !! Terminates visualiser + !! + !! Cleans up remnants of visualiser once it is no longer needed + !! + !! Result: + !! An empty visualiser object + !! + subroutine kill(self) + class(visualiser), intent(inout) :: self + + self % name ='' + self % geom => null() + call self % vizDict % kill() + + end subroutine kill + + + !! + !! Convert matIdx to a 24bit colour + !! + !! Special materials are associated with special colours: + !! OUTSIDE_MAT -> white (#ffffff) + !! VOID_MAT -> black (#000000) + !! UNDEF_MAT -> green (#00ff00) + !! + !! Args: + !! matIdx [in] -> Value of the material index + !! offset [in] -> Offset to be used in the hash function + !! + !! Result: + !! A 24-bit colour specifying the material + !! + elemental function materialColour(matIdx, offset) result(colour) + integer(shortInt), intent(in) :: matIdx + integer(shortInt) :: colour + integer(shortInt), intent(in) :: offset + + ! Since Knuth hash is cheap we can compute it anyway even if colour from + ! the map will end up being used + colour = mm_colourMap % getOrDefault(matIdx, knuthHash(matIdx + offset, 24)) + + end function materialColour + + !! + !! Convert uniqueID to 24bit colour + !! + !! An elemental wrapper over Knuth Hash + !! + !! We use a hash function to scatter colours across all available. + !! Knuth multiplicative hash is very good at scattering integer + !! sequences e.g. {1, 2, 3...}. Thus, it is ideal for a colourmap. + !! + !! Args: + !! uniqueID [in] -> Value of the uniqueID + !! + !! Result: + !! A 24-bit colour specifying the uniqueID + !! + elemental function uniqueIDColour(uniqueID) result(colour) + integer(shortInt), intent(in) :: uniqueID + integer(shortInt) :: colour + + colour = knuthHash(uniqueID, 24) + + end function uniqueIDColour + + +end module visualiser_class diff --git a/docs/Output.rst b/docs/Output.rst index 7349deb8d..9d3108e3b 100644 --- a/docs/Output.rst +++ b/docs/Output.rst @@ -38,8 +38,9 @@ Writing to an output file in SCONE is done through a series of calls to appropri character(nameLen) :: name type(outputFile) :: out - !! Initialise output by choosing the format - call out % init('asciiMatlab') + !! Initialise output by choosing the format and output file + !! If no file name is provided, the output will not be printed + call out % init('asciiMatlab', filename="./path/to/output/without_any_extension") !! Write an integer of 7 name = 'int' @@ -76,9 +77,10 @@ Writing to an output file in SCONE is done through a series of calls to appropri !! Same goes for blocks call out % endBlock() - !! When done, we can write output to a file - !! Extension will be determined by the format printer - call out % writeToFile("./path/to/output/without_any_extension") + !! We need to make sure that output file is properly finalised. We can skip + !! this line if we are not in a `program` block e.g. in `subroutine` or `function` + call out % finalise() + The need to write name to a `character(nameLen)` variable before calling the procedure is a quirk caused by the fact that the output file expects character of `nameLen` as a variable. If one were to diff --git a/docs/User Manual.rst b/docs/User Manual.rst index 65da84db2..7cddba610 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -275,9 +275,9 @@ neutronCEstd, to perform analog collision processing target nuclide movement. Target movement is sampled if target mass A < massThreshold. [Mn] * DBRCeMin (*optional*, default = 1.0e-08): minimum DBRC energy. [MeV] * DBRCeMax (*optional*, default = 2.0e-04): maximum DBRC energy. [MeV] - + Example: :: - + collisionOperator { neutronCE { type neutronCEstd; minEnergy 1.0e-12; maxEnergy 30.0; energyThreshold 200; massThreshold 2; DBRCeMin 1.0e-06; DBRCeMax 0.001; } } @@ -301,6 +301,7 @@ neutronCEimp, to perform implicit collision processing * minWgt (*optional*, default = 0.25): minimum particle weight for rouletting * maxWgt (*optional*, default = 1.25): maximum particle weight for splitting * avgWgt (*optional*, default = 0.5): weight of a particle on surviving rouletting +* maxSplit (*optional*, default = 1000): maximum number of splits allowed per particle * impAbs (*optional*, default = 0): 1 for true; 0 for false; enables implicit capture * impGen (*optional*, default = 1): 1 for true; 0 for false; enables implicit fission sites generation @@ -327,6 +328,19 @@ Example: :: collisionOperator { neutronMG { type neutronMGstd; } } +neutronMGimp +############ + +neutronMGimp, to perform implicit collision processing + +* maxSplit (*optional*, default = 1000): maximum number of splits allowed per particle +* weightWindows (*optional*, default = 0): 1 for true; 0 for false; enables the use of + weight windows + +Example: :: + + collisionOperator { neutronMG { type neutronMGimp; weightWindows 1; maxSplit 50; } } + Weight Windows -------------- @@ -679,6 +693,9 @@ bmp * what (*optional*, default = material): defines what is highlighted in the plot; options are ``material`` and ``uniqueID``, where ``uniqueID`` highlights unique cell IDs +* offset (*optional*, default = random) An integer (positive or negative) that + shifts the sequence of colours assigned to materials. Allows to change colours + from the default sequence in a parametric way. Example: :: @@ -736,7 +753,7 @@ Example: :: .. note:: If DBRC is applied, the 0K cross section ace files of the relevant nuclides must be included in the aceLibrary file. - + baseMgNeutronDatabase ##################### @@ -792,6 +809,9 @@ Other options are: * xsFile: needed for multi-group simulations. Must contain the path to the file where the multi-group cross sections are stored. +* rgb (*optional*): An array of three integers specifying the RGB colour e.g. ``(255 0 0)``. The + colour defined in this way will be used for visualisation of the material in the geometry plots. + Example 1: :: materials { @@ -802,6 +822,7 @@ Example 1: :: 8016.03 0.018535464; } } water { temp 273; + rgb (0 0 200); composition { 1001.03 0.0222222; 8016.03 0.00535; } @@ -1075,7 +1096,7 @@ Example: :: } } -.. note:: +.. note:: To calculate the average weight, one should divide weight moment 1 (weight1) by weight moment 0 (weight0). To calculate the variance of the weights, the tally results have to be post-processed as: var = weight2/weight0 - (weight1/weight0)^2