diff --git a/CMakeLists.txt b/CMakeLists.txt index a08c3b097..4ff09663d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,6 +77,9 @@ list(APPEND CMAKE_PREFIX_PATH $ENV{LAPACK_INSTALL}) find_package(LAPACK REQUIRED ) message(STATUS ${LAPACK_LIBRARIES}) +set(BLAS_LIBS /home/ajb343/BLAS/BLAS-3.10.0/Build) +find_package(BLAS REQUIRED) + # Dependencies for BUILD_TESTS if (BUILD_TESTS) # FIND PYTHON INTERPRETER diff --git a/CollisionOperator/CollisionProcessors/CMakeLists.txt b/CollisionOperator/CollisionProcessors/CMakeLists.txt index ee923a0ed..463788d8d 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 + ./IMCMGstd_class.f90) diff --git a/CollisionOperator/CollisionProcessors/IMCMGstd_class.f90 b/CollisionOperator/CollisionProcessors/IMCMGstd_class.f90 new file mode 100644 index 000000000..a86b2525e --- /dev/null +++ b/CollisionOperator/CollisionProcessors/IMCMGstd_class.f90 @@ -0,0 +1,220 @@ +module IMCMGstd_class + + use numPrecision + use endfConstants + 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_PHOTON + use particleDungeon_class, only : particleDungeon + + ! Abstract interface + use collisionProcessor_inter, only : collisionProcessor, collisionData ,init_super => init + + ! Nuclear Data Interface + use nuclearDataReg_mod, only : ndReg_getIMCMG => getIMCMG + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase + use mgIMCMaterial_inter, only : mgIMCMaterial, mgIMCMaterial_CptrCast + use reactionHandle_inter, only : reactionHandle + + implicit none + private + + !! + !! Standard (default) scalar collision processor for MG IMC + !! Determines type of collision as either absorption or effective scattering + !! + !! Settings: + !! NONE + !! + !! Sample dictionary input: + !! collProcName { + !! type IMCMGstd; + !! } + !! + type, public, extends(collisionProcessor) :: IMCMGstd + private + class(mgIMCDatabase), pointer, public :: xsData => null() + class(mgIMCMaterial), pointer, public :: mat => null() + contains + ! Initialisation procedure + procedure :: init + + ! Implementation of customisable procedures + procedure :: sampleCollision + procedure :: implicit + procedure :: elastic + procedure :: inelastic + procedure :: capture + procedure :: fission + procedure :: cutoffs + end type IMCMGstd + +contains + + !! + !! Initialise from dictionary + !! + subroutine init(self, dict) + class(IMCMGstd), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(100), parameter :: Here = 'init (IMCMGstd_class.f90)' + + ! Call superclass + call init_super(self, dict) + + end subroutine init + + !! + !! Samples collision + !! + !! Absorption with probability equal to fleck factor, otherwise + !! effective scattering + !! + !! Physical scattering is omitted as in reference paper "Four Decades of Implicit Monte Carlo" + !! (Allan B Wollaber) but may be included later if desired + !! + subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + real(defReal) :: r, fleck + character(100),parameter :: Here =' sampleCollision (IMCMGstd_class.f90)' + + ! Verify that particle is MG PHOTON + if( .not. p % isMG .or. p % type /= P_PHOTON) then + call fatalError(Here, 'Supports only MG PHOTON. Was given NEUTRON and/or CE '//printType(p % type)) + end if + + ! Verify and load nuclear data pointer + self % xsData => ndReg_getIMCMG() + if(.not.associated(self % xsData)) call fatalError(Here, "Failed to get active database for MG IMC") + + ! Get and verify material pointer + self % mat => mgIMCMaterial_CptrCast( self % xsData % getMaterial( p % matIdx())) + if(.not.associated(self % mat)) call fatalError(Here, "Failed to get MG IMC Material") + + r = p % pRNG % get() + + fleck = self % mat % getFleck() + + if( r < fleck ) then + ! Effective absoprtion + collDat % MT = macroCapture + else + ! Effective scattering + collDat % MT = macroAllScatter + end if + + end subroutine sampleCollision + + !! + !! Perform implicit treatment + !! + subroutine implicit(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + ! Do nothing. + + end subroutine implicit + + !! + !! Effective scattering - currently only elastic (constant energy-weight) + !! + subroutine elastic(self, p , collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + real(defReal) :: phi, mu + real(defReal), dimension(3) :: dir + character(100), parameter :: Here = 'elastic (IMCMGstd_class.f90)' + + ! Assign MT number + collDat % MT = macroAllScatter + + ! Sample Direction - chosen uniformly inside unit sphere + mu = 2 * p % pRNG % get() - 1 + phi = p % pRNG % get() * 2*pi + dir(1) = mu + dir(2) = sqrt(1-mu**2) * cos(phi) + dir(3) = sqrt(1-mu**2) * sin(phi) + + !p % coords % dir = dir + call p % rotate(mu, phi) + + end subroutine elastic + + !! + !! Inelastic scattering - Not currently supported + !! + subroutine inelastic(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + character(100),parameter :: Here = "inelastic (IMCMGstd_class.f90)" + + ! Do nothing. Should not be called + + call fatalError(Here, "Inelastic subroutine should not be called") + + end subroutine inelastic + + !! + !! Perform capture + !! + subroutine capture(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + p % isDead = .true. + + end subroutine capture + + !! + !! Perform fission + !! + subroutine fission(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + character(100), parameter :: Here = 'fission (IMCMGstd_class.f90)' + + ! Do nothing. Should not be called + + call fatalError(Here, "Fission subroutine should not be called") + + end subroutine fission + + !! + !! Apply cutoffs or post-collision implicit treatment + !! + subroutine cutoffs(self, p, collDat, thisCycle, nextCycle) + class(IMCMGstd), intent(inout) :: self + class(particle), intent(inout) :: p + type(collisionData), intent(inout) :: collDat + class(particleDungeon),intent(inout) :: thisCycle + class(particleDungeon),intent(inout) :: nextCycle + + ! Do nothing + + end subroutine cutoffs + +end module IMCMGstd_class diff --git a/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 b/CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90 index 3f190f3a4..ce45809e7 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 IMCMGstd_class, only : IMCMGstd implicit none private @@ -24,7 +25,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',& + 'IMCMGstd '] contains @@ -59,6 +61,10 @@ subroutine new_collisionProcessor(new,dict) allocate(neutronMGstd :: new) call new % init(dict) + case('IMCMGstd') + allocate(IMCMGstd :: new) + call new % init(dict) + !*** NEW COLLISION PROCESSOR TEMPLATE ***! !case('') ! allocate( :: new) diff --git a/CollisionOperator/collisionOperator_class.f90 b/CollisionOperator/collisionOperator_class.f90 index 07ee57a49..7c548a830 100644 --- a/CollisionOperator/collisionOperator_class.f90 +++ b/CollisionOperator/collisionOperator_class.f90 @@ -83,6 +83,11 @@ subroutine init(self, dict) self % lookupTable(P_MG, P_NEUTRON) = 2 end if + if(dict % isPresent('photonMG')) then + call new_collisionProcessor(self % physicsTable(3) % proc, dict % getDictPtr('photonMG')) + self % lookupTable(P_MG, P_PHOTON) = 3 + end if + end subroutine init !! diff --git a/DataStructures/dynArray_class.f90 b/DataStructures/dynArray_class.f90 index 09a648cf9..4913bc7e6 100644 --- a/DataStructures/dynArray_class.f90 +++ b/DataStructures/dynArray_class.f90 @@ -34,6 +34,7 @@ module dynArray_class procedure :: pop => pop_shortInt procedure :: empty => empty_shortInt procedure :: kill => kill_shortInt + procedure :: isPresent => isPresent_shortInt ! Private procedures @@ -248,4 +249,21 @@ pure subroutine kill_shortInt(self) end subroutine kill_shortInt + !! + !! Checks if item is present in array + !! + pure function isPresent_shortInt(self, item) result(isPresent) + class(dynIntArray), intent(in) :: self + integer(shortInt), intent(in) :: item + logical(defBool) :: isPresent + + isPresent = .false. + + ! Return false if array is empty + if (self % mySize == 0) return + + if (any(self % array(1:self % mySize) == item)) isPresent = .true. + + end function isPresent_shortInt + end module dynArray_class diff --git a/Geometry/CMakeLists.txt b/Geometry/CMakeLists.txt index 162a061c6..f459c0965 100644 --- a/Geometry/CMakeLists.txt +++ b/Geometry/CMakeLists.txt @@ -9,6 +9,7 @@ add_sources( ./csg_class.f90 ./geometry_inter.f90 ./geometryStd_class.f90 ./geometryReg_mod.f90 + ./discretiseGeom_class.f90 ) add_unit_tests( ./Tests/geomGraph_test.f90 diff --git a/Geometry/Cells/cellShelf_class.f90 b/Geometry/Cells/cellShelf_class.f90 index 58e29c926..2f556e3f9 100644 --- a/Geometry/Cells/cellShelf_class.f90 +++ b/Geometry/Cells/cellShelf_class.f90 @@ -47,13 +47,13 @@ module cellShelf_class !! Interface: !! init -> Initialise from a dictionary & surfaceShelf !! getPtr -> Get pointer to a cell given by its index - !! getIdx -> Return index of a cell fivent its ID + !! getIdx -> Return index of a cell given its ID !! getID -> Return cell ID given its index !! getFill -> Return content of the cell. If -ve it is universe ID. If +ve it is matIdx. !! getSize -> Return the number of cells (max cellIdx) !! kill -> Return to uninitialised state !! - !! NOTE: Becouse cells are stored as pointers, calling `kill` is crucial to prevent + !! NOTE: Because cells are stored as pointers, calling `kill` is crucial to prevent !! memory leaks. TODO: Add `final` procedure here ? !! type, public :: cellShelf diff --git a/Geometry/Universes/latUniverse_class.f90 b/Geometry/Universes/latUniverse_class.f90 index 798ebd9d8..a970a93c8 100644 --- a/Geometry/Universes/latUniverse_class.f90 +++ b/Geometry/Universes/latUniverse_class.f90 @@ -90,6 +90,7 @@ module latUniverse_class procedure :: distance procedure :: cross procedure :: cellOffset + procedure :: getSizeN end type latUniverse contains @@ -364,6 +365,17 @@ function cellOffset(self, coords) result (offset) end function cellOffset + !! + !! Return dimensions of lattice + !! + function getSizeN(self) result(sizeN) + class(latUniverse), intent(in) :: self + integer(shortInt), dimension(3) :: sizeN + + sizeN = self % sizeN + + end function getSizeN + !! !! Return to uninitialised state !! diff --git a/Geometry/discretiseGeom_class.f90 b/Geometry/discretiseGeom_class.f90 new file mode 100644 index 000000000..f017e17c7 --- /dev/null +++ b/Geometry/discretiseGeom_class.f90 @@ -0,0 +1,287 @@ +!! +!! Module to help with simplify the process of writing input files, specifically aimed at IMC but +!! may be useful for other applications if needed +!! +module discretiseGeom_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, numToChar + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + + ! Geometry + use geometry_inter, only : geometry + use geometryReg_mod, only : gr_geomPtr => geomPtr, & + gr_addGeom => addGeom, & + gr_geomIdx => geomIdx, & + gr_kill => kill + ! Nuclear Data + use materialMenu_mod, only : mm_matTemp => matTemp ,& + mm_matFile => matFile + use nuclearDataReg_mod, only : ndReg_init => init ,& + ndReg_activate => activate ,& + ndReg_kill => kill, & + ndReg_get => get + use nuclearDatabase_inter, only : nuclearDatabase + + ! Useful variables + class(geometry), pointer :: inputGeom + class(nuclearDatabase), pointer :: inputNucData + integer(shortInt), dimension(:), allocatable :: sizeN + +contains + + !! + !! Initialises geometry and nuclear database given in input file, and uses this as a baseline to + !! write new geometry and nuclear data dictionaries by splitting the input geometry into a + !! uniform grid. Makes new text files 'newGeom.txt' and 'newData.txt'. + !! + !! New geometry is a lattice universe, with each lattice cell containing the same material at its + !! centre in the initial input geometry, but with each cell initialising a new instance of this + !! material object, allowing for spatial temperature variation. + !! + !! For N lattice cells, will write to files N instances of: + !! In newGeom.txt => pj { id j; type pinUniverse; radii (0); fills (mj); } + !! In newData.txt => mj { temp ...; composition {} xsFile ./.../...; volume ...; } + !! (filled in with correct data for underlying material) + !! + !! In future, might be able to somehow modify the way that geometry and materials are built such + !! that separate instances of each mat are not needed + !! + !! Sample input: + !! dicretise { dimensions (10 1 1); } + !! + !! Args: + !! dict [in] -> Input file dictionary + !! newGeom [out] -> dictionary containing new geometry input + !! newData [out] -> dictionary containing new nuclear database input + !! + !! Errors: + !! invalid input dimensions + !! + subroutine discretise(dict, newGeom, newData) + class(dictionary), intent(in) :: dict + type(dictionary), intent(out) :: newGeom + type(dictionary), intent(out) :: newData + class(dictionary), pointer :: tempDict + character(nameLen) :: dataName, geomName + integer(shortInt) :: inputGeomIdx + character(100), parameter :: Here = 'discretise (discretiseGeom_class.f90)' + + ! Get discretisation settings + tempDict => dict % getDictPtr('discretise') + call tempDict % get(sizeN, 'dimensions') + if (any(sizeN <= 0)) call fatalError(Here, 'Invalid dimensions given') + + ! Build Nuclear Data using input + call ndReg_init(dict % getDictPtr("nuclearData")) + + ! Build geometry using input + tempDict => dict % getDictPtr('geometry') + geomName = 'inputGeom' + call gr_addGeom(geomName, tempDict) + inputGeomIdx = gr_geomIdx(geomName) + inputGeom => gr_geomPtr(inputGeomIdx) + + ! Activate input Nuclear Data + dataName = 'mg' + call ndReg_activate(P_PHOTON_MG, dataName, inputGeom % activeMats()) + inputNucData => ndReg_get(P_PHOTON_MG) + + ! Create files for materials and geometry + open(unit = 21, file = 'newGeom.txt') + open(unit = 22, file = 'newData.txt') + + ! Write required information to files + call writeToFiles(tempDict) + + close(21) + close(22) + + ! Kill input geom and nuclear database + call inputGeom % kill() + call inputNucData % kill() + call gr_kill() + call ndReg_kill() + + ! Create geometry and nuclear database from new files + call fileToDict(newData, './newData.txt') + call fileToDict(newGeom, './newGeom.txt') + + end subroutine discretise + + !! + !! Automatically write required information to files in the following format: + !! + !! newGeom.txt: + !! type geometryStd; + !! boundary (...); + !! graph {type shrunk;} + !! surfaces { outer {id 1; type box; origin (...); halfwidth (...);}} + !! cells {} + !! universes { + !! root {id rootID; type rootUniverse; border 1; fill u;} + !! lattice {id latID; type latUniverse; origin (0 0 0); pitch (...); shape (nx ny nz); padMat void; + !! map ( + !! 1 + !! 2 + !! . + !! . + !! . + !! nx*ny*nz); } + !! p1 {id 1; type pinUniverse; radii (0); fills (m1);} + !! p2 {id 2; type pinUniverse; radii (0); fills (m2);} + !! . + !! . + !! . + !! pN {id N; type pinUniverse; radii (0); fills (mN);} } + !! + !! newData.txt: + !! handles {mg {type baseMgIMCDatabase;}} + !! materials { + !! m1 {temp ...; composition {} xsFile ./...; volume ...;} + !! m2 {temp ...; composition {} xsFile ./...; volume ...;} + !! . + !! . + !! . + !! mN {temp ...; compositino {} xsFile ./...; volume ...;} } + !! + !! + !! Args: + !! dict [in] -> geometry input dictionary + !! + subroutine writeToFiles(dict) + class(dictionary), intent(in), pointer :: dict + character(nameLen) :: geomType, graphType + integer(shortInt), dimension(:), allocatable :: boundary + class(dictionary), pointer :: tempDict + real(defReal) :: volume + real(defReal), dimension(3) :: origin, halfwidth, pitch, corner, r + real(defReal), dimension(6) :: bounds + integer(shortInt) :: i, N, rootID, latID, matIdx, uniqueID + integer(shortInt), dimension(3) :: ijk + + ! Get bounds, pitch and lower corner + bounds = inputGeom % bounds() + do i = 1, 3 + pitch(i) = (bounds(i+3)-bounds(i)) / sizeN(i) + corner(i) = bounds(i) + end do + + ! Calculate volume of each cell and number of cells + volume = pitch(1) * pitch(2) * pitch(3) + N = sizeN(1) * sizeN(2) * sizeN(3) + + ! Get geometry settings from input file + call dict % get(geomType, 'type') + call dict % get(boundary, 'boundary') + tempDict => dict % getDictPtr('graph') + call tempDict % get(graphType, 'type') + + ! Write to new file + write (21, '(8A)') 'type '//trim(geomType)//';'//new_line('A')//& + &'boundary ('//numToChar(boundary)//');'//new_line('A')//& + &'graph {type '//trim(graphType)//';}' + + ! Construct outer boundary surface based on input geometry bounds + ! Will likely cause problems if input outer surface is not a box + do i=1, 3 + origin(i) = (bounds(i+3) + bounds(i)) / 2 + halfwidth(i) = (bounds(i+3) - bounds(i)) / 2 + end do + write (21, '(8A)') 'surfaces { outer {id 1; type box; origin ('//numToChar(origin)//'); & + &halfwidth ('//numToChar(halfwidth)//');}}' + + ! Empty cell dictionary and construct root universe + rootID = N + 1 + latID = N + 2 + write (21, '(8A)') 'cells {}'//new_line('A')//& + &'universes {'//new_line('A')//& + &'root {id '//numToChar(rootID)//'; type rootUniverse; border 1; & + &fill u<'//numToChar(latID)//'>;}' + + ! Write lattice input + write(21, '(8A)') 'lattice {id '//numToChar(latID)//'; type latUniverse; origin (0 0 0); pitch ('//& + &numToChar(pitch)//'); shape ('//numToChar(sizeN)//'); padMat void;'& + &//new_line('A')//'map (' + do i = 1, N + write(21, '(8A)') numToChar(i) + end do + + write(21, '(8A)') '); }'//new_line('A') + + + ! Write material settings - TODO obtain this from input file automatically? + write (22, '(8A)') 'handles {mg {type baseMgIMCDatabase;}}'//new_line('A')//'materials {' + + + ! Write pin universes and materials + do i = 1, N + + ! Get location in centre of lattice cell + ijk = get_ijk(i, sizeN) + r = corner + pitch * (ijk - HALF) + + ! Get matIdx from input geometry + call inputGeom % whatIsAt(matIdx, uniqueId, r) + + ! Pin universe for void and outside mat + if (matIdx == VOID_MAT) then + write(21, '(8A)') 'p'//numToChar(i)//' {id '//numToChar(i)//'; type pinUniverse; & + &radii (0); fills (void);}' + + else if (matIdx == OUTSIDE_MAT) then + write(21, '(8A)') 'p'//numToChar(i)//' {id '//numToChar(i)//'; type pinUniverse; & + &radii (0); fills (outside);}' + + ! User-defined materials + else + ! Pin universe + write(21, '(8A)') 'p'//numToChar(i)//' {id '//numToChar(i)//'; type pinUniverse; & + &radii (0); fills (m'//numToChar(i)//');}' + + ! Material + write(22, '(8A)') 'm'//numToChar(i)//' {temp '//numToChar(mm_matTemp(matIdx))//'; & + &composition {} xsFile '//trim(mm_matFile(matIdx))//'; volume '//numToChar(volume)//';}' + + end if + end do + + ! Write closing brackets for dictionaries + write (21, '(8A)') '}' + write (22, '(8A)') '}' + + end subroutine writeToFiles + + !! + !! Generate ijk from localID and shape + !! + !! Args: + !! localID [in] -> Local id of the cell between 1 and product(sizeN) + !! sizeN [in] -> Number of cells in each cardinal direction x, y & z + !! + !! Result: + !! Array ijk which has integer position in each cardinal direction + !! + pure function get_ijk(localID, sizeN) result(ijk) + integer(shortInt), intent(in) :: localID + integer(shortInt), dimension(3), intent(in) :: sizeN + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: temp, base + + temp = localID - 1 + + base = temp / sizeN(1) + ijk(1) = temp - sizeN(1) * base + 1 + + temp = base + base = temp / sizeN(2) + ijk(2) = temp - sizeN(2) * base + 1 + + ijk(3) = base + 1 + + end function get_ijk + + +end module discretiseGeom_class diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index 09e59eee0..06ccf5c54 100644 --- a/Geometry/geometryStd_class.f90 +++ b/Geometry/geometryStd_class.f90 @@ -9,6 +9,7 @@ module geometryStd_class use geometry_inter, only : geometry, distCache use csg_class, only : csg use universe_inter, only : universe + use latUniverse_class, only : latUniverse use surface_inter, only : surface ! Nuclear Data @@ -57,6 +58,7 @@ module geometryStd_class procedure :: moveGlobal procedure :: teleport procedure :: activeMats + procedure :: latSizeN ! Private procedures procedure, private :: diveToMat @@ -572,5 +574,42 @@ subroutine closestDist_cache(self, dist, surfIdx, lvl, coords, cache) end do end subroutine closestDist_cache + !! + !! Return dimensions of latUniverse + !! + !! fatalError if no latUniverse found, if there are multiple then it will return dimensions + !! of the first one found, which may not be what is wanted + !! + function latSizeN(self) result(sizeN) + class(geometryStd), intent(in) :: self + integer(shortInt), dimension(3) :: sizeN + integer(shortInt) :: i + class(universe), pointer :: uni + class(latUniverse), pointer :: latUni + character(100), parameter :: Here = 'latSizeN (geometryStd_class.f90)' + + ! Search for latUniverse + do i=1, self % geom % unis % getSize() + + uni => self % geom % unis % getPtr(i) + + select type(uni) + class is(latUniverse) + latUni => uni + exit + + class default + latUni => null() + + end select + end do + + if (.not. associated(latUni)) call fatalError(Here, 'Lattice universe not found') + + ! Find lattice dimensions + sizeN = latUni % getSizeN() + + end function latSizeN + end module geometryStd_class diff --git a/Geometry/geometry_inter.f90 b/Geometry/geometry_inter.f90 index 973007549..eb4370248 100644 --- a/Geometry/geometry_inter.f90 +++ b/Geometry/geometry_inter.f90 @@ -43,6 +43,7 @@ module geometry_inter procedure(moveGlobal), deferred :: moveGlobal procedure(teleport), deferred :: teleport procedure(activeMats), deferred :: activeMats + procedure(latSizeN), deferred :: latSizeN ! Common procedures procedure :: slicePlot @@ -262,6 +263,18 @@ function activeMats(self) result(matList) integer(shortInt), dimension(:), allocatable :: matList end function activeMats + !! + !! Return dimensions of latUniverse + !! + !! fatalError if no latUniverse found, if there are multiple then it will return dimensions + !! of the first one found, which may not be what is wanted + !! + function latSizeN(self) result(sizeN) + import geometry, shortInt + class(geometry), intent(in) :: self + integer(shortInt), dimension(3) :: sizeN + end function latSizeN + end interface contains diff --git a/InputFiles/IMC/DataFiles/hohlraumData b/InputFiles/IMC/DataFiles/hohlraumData new file mode 100644 index 000000000..1cdf6c285 --- /dev/null +++ b/InputFiles/IMC/DataFiles/hohlraumData @@ -0,0 +1,22 @@ +// This XSS are for a URR-2-1-IN/SL Benchamrk Problem +// Source: A. Sood, R. A. Forster, and D. K. Parsons, +// ‘Analytical Benchmark Test Set For Criticality Code Verification’ +// + +numberOfGroups 1; + +capture ( + 100 + -3 +); + +scatter ( + 0 + 0 +); + +cv ( + 0.3 //3.00026706e15 + 0 +); + diff --git a/InputFiles/IMC/DataFiles/marshakData b/InputFiles/IMC/DataFiles/marshakData new file mode 100644 index 000000000..389b9eaa0 --- /dev/null +++ b/InputFiles/IMC/DataFiles/marshakData @@ -0,0 +1,18 @@ + +numberOfGroups 1; + +capture ( + 10 + -3 +); + +scatter ( + 0 + 0 +); + +cv ( + 7.14 + 0 +); + diff --git a/InputFiles/IMC/DataFiles/sampleData b/InputFiles/IMC/DataFiles/sampleData new file mode 100644 index 000000000..248608a6e --- /dev/null +++ b/InputFiles/IMC/DataFiles/sampleData @@ -0,0 +1,44 @@ +// +// Sample material data file for IMC calculations +// + + +numberOfGroups 1; + +alpha = 1; + // Optional setting of alpha used in calculation of fleck factor. If not given, alpha is set at 1. + + + // Set polynomial temperature-dependent opacities for the material. + // Currently have only considered the grey case, if using a frequency dependent opacity + // then this would need to be changed to a more complex input. + // Input should be a 1D array of coefficients followed by exponents, with any polynomial + // length allowed + // e.g. Here, sigmaA = 1 + 2T + +capture ( // Absorption opacity + 1 2 // Coefficients + 0 1 // Exponents +); + +scatter ( // Scattering opacity + 0 + 0 +); + + + // Set temperature-dependent specific heat capacity of the material. + // Same format as above. + // Currently cannot have an exponent of exactly -1, as cv is integrated simply by adding 1 to + // exponents an have not yet allowed T^(-1) to integrate to ln(T). + // After integration, solved by Newton-Raphson solver. Some choices of cv may not converge, + // and some will give negative energies and temperatures. Unsure if this is due to some + // numerical oversight in the way the calculation is done or if these are just unphysical + // choices. + // e.g. Here, cv = 4T^3 - 2T + T^(-0.5) + +cv ( + 4 -2 1 // Coefficients + 3 1 -0.5 // Exponents +); + diff --git a/InputFiles/IMC/SimpleCases/3region b/InputFiles/IMC/SimpleCases/3region new file mode 100644 index 000000000..4a17874b1 --- /dev/null +++ b/InputFiles/IMC/SimpleCases/3region @@ -0,0 +1,79 @@ + +type IMCPhysicsPackage; + +pop 5000; +limit 20000; +steps 50; +timeStepSize 0.1; +printUpdates 3; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +tally { + } + +geometry { + type geometryStd; + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + sep1 { id 1; type xPlane; x0 -0.5; } + sep2 { id 2; type xPlane; x0 0.5; } + outer { id 3; type box; origin ( 0.0 0.0 0.0); halfwidth ( 1.5 0.5 0.5); } + } + cells + { + cell1 { id 1; type simpleCell; surfaces ( -1); filltype mat; material mat1; } + cell2 { id 2; type simpleCell; surfaces (1 -2); filltype mat; material mat2; } + cell3 { id 3; type simpleCell; surfaces ( 2); filltype mat; material mat3; } + } + universes + { + root { id 1; type rootUniverse; border 3; fill u<2>; } + cell { id 2; type cellUniverse; cells ( 1 2 3); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { + temp 1; + composition {} + xsFile ./dataFiles/imcData; + volume 1; + } + mat2 { + temp 0.01; + composition {} + xsFile ./dataFiles/imcData; + volume 1; + } + mat3 { + temp 0.01; + composition {} + xsFile ./dataFiles/imcData; + volume 1; + } + + +} + +} + + + diff --git a/InputFiles/IMC/SimpleCases/dataFiles/imcData b/InputFiles/IMC/SimpleCases/dataFiles/imcData new file mode 100644 index 000000000..3e5ec6ed5 --- /dev/null +++ b/InputFiles/IMC/SimpleCases/dataFiles/imcData @@ -0,0 +1,18 @@ + +numberOfGroups 1; + +capture ( + 1 + 0 +); + +scatter ( + 0 + 0 +); + +cv ( + 4 + 3 +); + diff --git a/InputFiles/IMC/SimpleCases/dataFiles/imcData2 b/InputFiles/IMC/SimpleCases/dataFiles/imcData2 new file mode 100644 index 000000000..0e7f39427 --- /dev/null +++ b/InputFiles/IMC/SimpleCases/dataFiles/imcData2 @@ -0,0 +1,18 @@ + +numberOfGroups 1; + +capture ( + 1.0 + 0.0 +); + +scatter ( + 0.0 + 0.0 +); + +cv ( + 4.0 3.0 + 3.0 2.0 +); + diff --git a/InputFiles/IMC/SimpleCases/infiniteRegion b/InputFiles/IMC/SimpleCases/infiniteRegion new file mode 100644 index 000000000..5ec1dde1c --- /dev/null +++ b/InputFiles/IMC/SimpleCases/infiniteRegion @@ -0,0 +1,65 @@ + +type IMCPhysicsPackage; + +pop 5000; +limit 20000; +steps 50; +timeStepSize 0.01; +printUpdates 1; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +tally { + } + +geometry { + type geometryStd; + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + squareBound { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 0.5 0.5 0.5); } + } + cells {} + universes + { + + root + { + id 1; + type rootUniverse; + border 1; + fill mat; + } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat { + temp 1; + composition {} + xsFile ./dataFiles/imcData; + volume 1; + } + +} + +} + + + diff --git a/InputFiles/IMC/SimpleCases/sphereInCube b/InputFiles/IMC/SimpleCases/sphereInCube new file mode 100644 index 000000000..556e53737 --- /dev/null +++ b/InputFiles/IMC/SimpleCases/sphereInCube @@ -0,0 +1,69 @@ + +type IMCPhysicsPackage; + +pop 100; +limit 2000; +steps 500; +timeStepSize 0.1; +printUpdates 2; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +tally { + } + +geometry { + type geometryStd; + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + inner { id 1; type sphere; origin ( 0.0 0.0 0.0); radius 1; } + outer { id 2; type box; origin ( 0.0 0.0 0.0); halfwidth ( 1 1 1); } + } + cells + { + inner_cell { id 1; type simpleCell; surfaces (-1); filltype mat; material mat1; } + outer_cell { id 2; type simpleCell; surfaces ( 1); filltype mat; material mat2; } + } + universes + { + root { id 1; type rootUniverse; border 2; fill u<2>; } + cell { id 2; type cellUniverse; cells ( 1 2 ); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { + temp 1; + composition {} + xsFile ./dataFiles/imcData; + volume 4.18879; + } + mat2 { + temp 5; + composition {} + xsFile ./dataFiles/imcData; + volume 3.81121; } + +} + +} + + + diff --git a/InputFiles/IMC/SimpleCases/touchingCubes b/InputFiles/IMC/SimpleCases/touchingCubes new file mode 100644 index 000000000..3cfbfa1fb --- /dev/null +++ b/InputFiles/IMC/SimpleCases/touchingCubes @@ -0,0 +1,70 @@ + +type IMCPhysicsPackage; + +pop 5000; +limit 20000; +steps 50; +timeStepSize 1; +printUpdates 2; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +tally { + } + +geometry { + type geometryStd; + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + sep { id 1; type xPlane; x0 0.0; } + outer { id 2; type box; origin ( 0.0 0.0 0.0); halfwidth ( 1 0.5 0.5); } + } + cells + { + cell1 { id 1; type simpleCell; surfaces (-1); filltype mat; material mat1; } + cell2 { id 2; type simpleCell; surfaces ( 1); filltype mat; material mat2; } + } + universes + { + root { id 1; type rootUniverse; border 2; fill u<2>; } + cell { id 2; type cellUniverse; cells ( 1 2 ); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { + temp 1; + composition {} + xsFile ./dataFiles/imcData; + volume 1; + } + mat2 { + temp 0; + composition {} + xsFile ./dataFiles/imcData2; + volume 1; + } + +} + +} + + + diff --git a/InputFiles/IMC/hohlraum b/InputFiles/IMC/hohlraum new file mode 100644 index 000000000..2be0bed98 --- /dev/null +++ b/InputFiles/IMC/hohlraum @@ -0,0 +1,100 @@ + +type IMCPhysicsPackage; + +pop 2000000; +limit 20000000; +steps 100; +timeStepSize 0.00000000001; +printUpdates 0; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + method ST; + cutoff 0.9; + grid { dimensions (20 20 1); searchN (10 10 1); } + } + +matSource { type imcSource; method fast; } + +source { + type bbSurfaceSource; + r (-0.5 -0.5 -0.5 0.5 -0.5 0.5); + temp 1; + N 20000; +} + +viz { vizDict { type vtk; corner (-0.5 -0.5 -0.5); width (1 1 1); vox (20 20 1); } } + +tally { + } + +discretise { dimensions (200 200 1); } + +geometry { + type geometryStd; + boundary (0 0 0 0 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 0.5 0.5 0.5); } + xPlane1 { id 2; type xPlane; x0 -0.45; } + xPlane2 { id 3; type xPlane; x0 -0.25; } + xPlane3 { id 4; type xPlane; x0 0.25; } + xPlane4 { id 5; type xPlane; x0 0.45; } + yPlane1 { id 6; type yPlane; y0 -0.45; } + yPlane2 { id 7; type yPlane; y0 -0.25; } + yPlane3 { id 8; type yPlane; y0 0.25; } + yPlane4 { id 9; type yPlane; y0 0.45; } + } + + cells { + // Top mat + cell1 { id 1; type simpleCell; surfaces (9); filltype mat; material mat1; } + // Bottom mat + cell2 { id 2; type simpleCell; surfaces (-6); filltype mat; material mat1; } + // Right side mat + cell3 { id 3; type simpleCell; surfaces (5 6 -9); filltype mat; material mat1; } + // Left side mat + cell4 { id 4; type simpleCell; surfaces (-2 7 -8); filltype mat; material mat1; } + // Centre mat + cell5 { id 5; type simpleCell; surfaces (3 -4 7 -8); filltype mat; material mat1; } + // Top void + cell6 { id 6; type simpleCell; surfaces (8 -9 -5); filltype mat; material void; } + // Bottom void + cell7 { id 7; type simpleCell; surfaces (6 -7 -5); filltype mat; material void; } + // Right void + cell8 { id 8; type simpleCell; surfaces (7 -8 4 -5); filltype mat; material void; } + // Left void + cell9 { id 9; type simpleCell; surfaces (7 -8 2 -3); filltype mat; material void; } + } + + universes + { + root { id 1; type rootUniverse; border 1; fill u<2>; } + cells { id 2; type cellUniverse; cells (1 2 3 4 5 6 7 8 9); } + } + +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.000000001; composition {} xsFile ./DataFiles/hohlraumData; } + + } + +} + + + diff --git a/InputFiles/IMC/marshakWave b/InputFiles/IMC/marshakWave new file mode 100644 index 000000000..c3d87c7ae --- /dev/null +++ b/InputFiles/IMC/marshakWave @@ -0,0 +1,49 @@ +// Marshak wave IMC benchmark +// Requires lightSpeed and radiation constant set to ONE in universalVariables + +type IMCPhysicsPackage; + +pop 1000000; +limit 2500000; +steps 10; +timeStepSize 0.5; +printUpdates 20; + +collisionOperator { photonMG {type IMCMGstd;} } + +transportOperator { + type transportOperatorTimeHT; cutoff 0.8; + } + +tally {} + +// Material photon source - optional but will default to rejection sampling if method not specified +matSource { type imcSource; method fast; } + +// Black body surface source +source { type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 10000; } + +// Spatial discretisation +discretise { dimensions (500 1 1); } + +// Overlaid grid for hybrid tracking +grid { dimensions (50 1 1); searchN (10 1 1); } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces { outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); }} + cells {} + universes { root { id 100; type rootUniverse; border 1; fill mat1; }} +} + +nuclearData { + handles { mg { type baseMgIMCDatabase; }} + + materials { mat1 { temp 0.01; composition {} xsFile ./DataFiles/marshakData; }} +} + + + diff --git a/InputFiles/IMC/oldInputs/dataFiles/imcData b/InputFiles/IMC/oldInputs/dataFiles/imcData new file mode 100644 index 000000000..389b9eaa0 --- /dev/null +++ b/InputFiles/IMC/oldInputs/dataFiles/imcData @@ -0,0 +1,18 @@ + +numberOfGroups 1; + +capture ( + 10 + -3 +); + +scatter ( + 0 + 0 +); + +cv ( + 7.14 + 0 +); + diff --git a/InputFiles/IMC/oldInputs/marshakWave128 b/InputFiles/IMC/oldInputs/marshakWave128 new file mode 100644 index 000000000..7f26f2a56 --- /dev/null +++ b/InputFiles/IMC/oldInputs/marshakWave128 @@ -0,0 +1,348 @@ +// Marshak wave simulation using 128 equal spatial regions +// *** lightSpeed and radiationConstant must both be changed to ONE in universalVariables *** + +type IMCPhysicsPackage; + +pop 10000; +limit 200000; +steps 1000; +timeStepSize 0.05; +printUpdates 4; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; cutoff 0.9; + } + +source { + type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 1000; +} + +tally { + } + +grid { dimensions (25 1 1); searchN (1000 1 1); } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); } + } + + cells + { + } + universes + { + root { id 1000; type rootUniverse; border 1; fill u<2000>; } + + lat { id 2000; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (0.03125 1.0 1.0); + shape (128 1 1); + padMat void; + map ( 1 2 3 4 5 6 7 8 9 10 + 11 12 13 14 15 16 17 18 19 20 + 21 22 23 24 25 26 27 28 29 30 + 31 32 33 34 35 36 37 38 39 40 + 41 42 43 44 45 46 47 48 49 50 + 51 52 53 54 55 56 57 58 59 60 + 61 62 63 64 65 66 67 68 69 70 + 71 72 73 74 75 76 77 78 79 80 + 81 82 83 84 85 86 87 88 89 90 + 91 92 93 94 95 96 97 98 99 100 + 101 102 103 104 105 106 107 108 109 110 + 111 112 113 114 115 116 117 118 119 120 + 121 122 123 124 125 126 127 128); + } + + zone1 { id 1; type pinUniverse; radii (0.0); fills (mat1); } + zone2 { id 2; type pinUniverse; radii (0.0); fills (mat2); } + zone3 { id 3; type pinUniverse; radii (0.0); fills (mat3); } + zone4 { id 4; type pinUniverse; radii (0.0); fills (mat4); } + zone5 { id 5; type pinUniverse; radii (0.0); fills (mat5); } + zone6 { id 6; type pinUniverse; radii (0.0); fills (mat6); } + zone7 { id 7; type pinUniverse; radii (0.0); fills (mat7); } + zone8 { id 8; type pinUniverse; radii (0.0); fills (mat8); } + zone9 { id 9; type pinUniverse; radii (0.0); fills (mat9); } + zone10 { id 10; type pinUniverse; radii (0.0); fills (mat10); } + zone11 { id 11; type pinUniverse; radii (0.0); fills (mat11); } + zone12 { id 12; type pinUniverse; radii (0.0); fills (mat12); } + zone13 { id 13; type pinUniverse; radii (0.0); fills (mat13); } + zone14 { id 14; type pinUniverse; radii (0.0); fills (mat14); } + zone15 { id 15; type pinUniverse; radii (0.0); fills (mat15); } + zone16 { id 16; type pinUniverse; radii (0.0); fills (mat16); } + zone17 { id 17; type pinUniverse; radii (0.0); fills (mat17); } + zone18 { id 18; type pinUniverse; radii (0.0); fills (mat18); } + zone19 { id 19; type pinUniverse; radii (0.0); fills (mat19); } + zone20 { id 20; type pinUniverse; radii (0.0); fills (mat20); } + zone21 { id 21; type pinUniverse; radii (0.0); fills (mat21); } + zone22 { id 22; type pinUniverse; radii (0.0); fills (mat22); } + zone23 { id 23; type pinUniverse; radii (0.0); fills (mat23); } + zone24 { id 24; type pinUniverse; radii (0.0); fills (mat24); } + zone25 { id 25; type pinUniverse; radii (0.0); fills (mat25); } + zone26 { id 26; type pinUniverse; radii (0.0); fills (mat26); } + zone27 { id 27; type pinUniverse; radii (0.0); fills (mat27); } + zone28 { id 28; type pinUniverse; radii (0.0); fills (mat28); } + zone29 { id 29; type pinUniverse; radii (0.0); fills (mat29); } + zone30 { id 30; type pinUniverse; radii (0.0); fills (mat30); } + zone31 { id 31; type pinUniverse; radii (0.0); fills (mat31); } + zone32 { id 32; type pinUniverse; radii (0.0); fills (mat32); } + + zone33 { id 33; type pinUniverse; radii (0.0); fills (mat33); } + zone34 { id 34; type pinUniverse; radii (0.0); fills (mat34); } + zone35 { id 35; type pinUniverse; radii (0.0); fills (mat35); } + zone36 { id 36; type pinUniverse; radii (0.0); fills (mat36); } + zone37 { id 37; type pinUniverse; radii (0.0); fills (mat37); } + zone38 { id 38; type pinUniverse; radii (0.0); fills (mat38); } + zone39 { id 39; type pinUniverse; radii (0.0); fills (mat39); } + zone40 { id 40; type pinUniverse; radii (0.0); fills (mat40); } + zone41 { id 41; type pinUniverse; radii (0.0); fills (mat41); } + zone42 { id 42; type pinUniverse; radii (0.0); fills (mat42); } + zone43 { id 43; type pinUniverse; radii (0.0); fills (mat43); } + zone44 { id 44; type pinUniverse; radii (0.0); fills (mat44); } + zone45 { id 45; type pinUniverse; radii (0.0); fills (mat45); } + zone46 { id 46; type pinUniverse; radii (0.0); fills (mat46); } + zone47 { id 47; type pinUniverse; radii (0.0); fills (mat47); } + zone48 { id 48; type pinUniverse; radii (0.0); fills (mat48); } + zone49 { id 49; type pinUniverse; radii (0.0); fills (mat49); } + zone50 { id 50; type pinUniverse; radii (0.0); fills (mat50); } + zone51 { id 51; type pinUniverse; radii (0.0); fills (mat51); } + zone52 { id 52; type pinUniverse; radii (0.0); fills (mat52); } + zone53 { id 53; type pinUniverse; radii (0.0); fills (mat53); } + zone54 { id 54; type pinUniverse; radii (0.0); fills (mat54); } + zone55 { id 55; type pinUniverse; radii (0.0); fills (mat55); } + zone56 { id 56; type pinUniverse; radii (0.0); fills (mat56); } + zone57 { id 57; type pinUniverse; radii (0.0); fills (mat57); } + zone58 { id 58; type pinUniverse; radii (0.0); fills (mat58); } + zone59 { id 59; type pinUniverse; radii (0.0); fills (mat59); } + zone60 { id 60; type pinUniverse; radii (0.0); fills (mat60); } + zone61 { id 61; type pinUniverse; radii (0.0); fills (mat61); } + zone62 { id 62; type pinUniverse; radii (0.0); fills (mat62); } + zone63 { id 63; type pinUniverse; radii (0.0); fills (mat63); } + zone64 { id 64; type pinUniverse; radii (0.0); fills (mat64); } + + zone65 { id 65; type pinUniverse; radii (0.0); fills (mat65); } + zone66 { id 66; type pinUniverse; radii (0.0); fills (mat66); } + zone67 { id 67; type pinUniverse; radii (0.0); fills (mat67); } + zone68 { id 68; type pinUniverse; radii (0.0); fills (mat68); } + zone69 { id 69; type pinUniverse; radii (0.0); fills (mat69); } + zone70 { id 70; type pinUniverse; radii (0.0); fills (mat70); } + zone71 { id 71; type pinUniverse; radii (0.0); fills (mat71); } + zone72 { id 72; type pinUniverse; radii (0.0); fills (mat72); } + zone73 { id 73; type pinUniverse; radii (0.0); fills (mat73); } + zone74 { id 74; type pinUniverse; radii (0.0); fills (mat74); } + zone75 { id 75; type pinUniverse; radii (0.0); fills (mat75); } + zone76 { id 76; type pinUniverse; radii (0.0); fills (mat76); } + zone77 { id 77; type pinUniverse; radii (0.0); fills (mat77); } + zone78 { id 78; type pinUniverse; radii (0.0); fills (mat78); } + zone79 { id 79; type pinUniverse; radii (0.0); fills (mat79); } + zone80 { id 80; type pinUniverse; radii (0.0); fills (mat80); } + zone81 { id 81; type pinUniverse; radii (0.0); fills (mat81); } + zone82 { id 82; type pinUniverse; radii (0.0); fills (mat82); } + zone83 { id 83; type pinUniverse; radii (0.0); fills (mat83); } + zone84 { id 84; type pinUniverse; radii (0.0); fills (mat84); } + zone85 { id 85; type pinUniverse; radii (0.0); fills (mat85); } + zone86 { id 86; type pinUniverse; radii (0.0); fills (mat86); } + zone87 { id 87; type pinUniverse; radii (0.0); fills (mat87); } + zone88 { id 88; type pinUniverse; radii (0.0); fills (mat88); } + zone89 { id 89; type pinUniverse; radii (0.0); fills (mat89); } + zone90 { id 90; type pinUniverse; radii (0.0); fills (mat90); } + zone91 { id 91; type pinUniverse; radii (0.0); fills (mat91); } + zone92 { id 92; type pinUniverse; radii (0.0); fills (mat92); } + zone93 { id 93; type pinUniverse; radii (0.0); fills (mat93); } + zone94 { id 94; type pinUniverse; radii (0.0); fills (mat94); } + zone95 { id 95; type pinUniverse; radii (0.0); fills (mat95); } + zone96 { id 96; type pinUniverse; radii (0.0); fills (mat96); } + + zone97 { id 97; type pinUniverse; radii (0.0); fills (mat97); } + zone98 { id 98; type pinUniverse; radii (0.0); fills (mat98); } + zone99 { id 99; type pinUniverse; radii (0.0); fills (mat99); } + zone100 { id 100; type pinUniverse; radii (0.0); fills (mat100); } + zone101 { id 101; type pinUniverse; radii (0.0); fills (mat101); } + zone102 { id 102; type pinUniverse; radii (0.0); fills (mat102); } + zone103 { id 103; type pinUniverse; radii (0.0); fills (mat103); } + zone104 { id 104; type pinUniverse; radii (0.0); fills (mat104); } + zone105 { id 105; type pinUniverse; radii (0.0); fills (mat105); } + zone106 { id 106; type pinUniverse; radii (0.0); fills (mat106); } + zone107 { id 107; type pinUniverse; radii (0.0); fills (mat107); } + zone108 { id 108; type pinUniverse; radii (0.0); fills (mat108); } + zone109 { id 109; type pinUniverse; radii (0.0); fills (mat109); } + zone110 { id 110; type pinUniverse; radii (0.0); fills (mat110); } + zone111 { id 111; type pinUniverse; radii (0.0); fills (mat111); } + zone112 { id 112; type pinUniverse; radii (0.0); fills (mat112); } + zone113 { id 113; type pinUniverse; radii (0.0); fills (mat113); } + zone114 { id 114; type pinUniverse; radii (0.0); fills (mat114); } + zone115 { id 115; type pinUniverse; radii (0.0); fills (mat115); } + zone116 { id 116; type pinUniverse; radii (0.0); fills (mat116); } + zone117 { id 117; type pinUniverse; radii (0.0); fills (mat117); } + zone118 { id 118; type pinUniverse; radii (0.0); fills (mat118); } + zone119 { id 119; type pinUniverse; radii (0.0); fills (mat119); } + zone120 { id 120; type pinUniverse; radii (0.0); fills (mat120); } + zone121 { id 121; type pinUniverse; radii (0.0); fills (mat121); } + zone122 { id 122; type pinUniverse; radii (0.0); fills (mat122); } + zone123 { id 123; type pinUniverse; radii (0.0); fills (mat123); } + zone124 { id 124; type pinUniverse; radii (0.0); fills (mat124); } + zone125 { id 125; type pinUniverse; radii (0.0); fills (mat125); } + zone126 { id 126; type pinUniverse; radii (0.0); fills (mat126); } + zone127 { id 127; type pinUniverse; radii (0.0); fills (mat127); } + zone128 { id 128; type pinUniverse; radii (0.0); fills (mat128); } + + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat2 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat3 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat4 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat5 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat6 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat7 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat8 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat9 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat10 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat11 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat12 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat13 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat14 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat15 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat16 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat17 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat18 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat19 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat20 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat21 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat22 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat23 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat24 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat25 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat26 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat27 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat28 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat29 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat30 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat31 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat32 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + + mat33 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat34 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat35 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat36 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat37 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat38 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat39 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat40 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat41 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat42 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat43 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat44 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat45 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat46 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat47 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat48 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat49 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat50 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat51 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat52 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat53 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat54 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat55 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat56 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat57 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat58 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat59 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat60 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat61 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat62 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat63 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat64 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + + mat65 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat66 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat67 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat68 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat69 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat70 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat71 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat72 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat73 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat74 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat75 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat76 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat77 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat78 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat79 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat80 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat81 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat82 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat83 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat84 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat85 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat86 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat87 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat88 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat89 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat90 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat91 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat92 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat93 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat94 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat95 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat96 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + + mat97 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat98 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat99 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat100 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat101 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat102 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat103 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat104 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat105 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat106 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat107 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat108 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat109 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat110 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat111 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat112 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat113 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat114 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat115 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat116 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat117 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat118 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat119 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat120 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat121 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat122 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat123 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat124 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat125 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat126 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat127 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + mat128 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.03125; } + + } + +} + + + diff --git a/InputFiles/IMC/oldInputs/marshakWave16 b/InputFiles/IMC/oldInputs/marshakWave16 new file mode 100644 index 000000000..6d4abf262 --- /dev/null +++ b/InputFiles/IMC/oldInputs/marshakWave16 @@ -0,0 +1,104 @@ +// Marshak wave simulation using 16 equal spatial regions +// *** lightSpeed and radiationConstant must both be changed to ONE in universalVariables *** + +type IMCPhysicsPackage; + +pop 16000; +limit 320000; +steps 1000; +timeStepSize 0.5; +printUpdates 4; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; cutoff 0.0; + } + +source { + type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 1000; +} + +tally { + } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); } + } + + cells + { + } + universes + { + root { id 100; type rootUniverse; border 1; fill u<200>; } + + lat { id 200; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (0.25 1.0 1.0); + shape (16 1 1); + padMat void; + map ( 1 2 3 4 5 6 7 8 9 10 + 11 12 13 14 15 16); + } + + zone1 { id 1; type pinUniverse; radii (0.0); fills (mat1); } + zone2 { id 2; type pinUniverse; radii (0.0); fills (mat2); } + zone3 { id 3; type pinUniverse; radii (0.0); fills (mat3); } + zone4 { id 4; type pinUniverse; radii (0.0); fills (mat4); } + zone5 { id 5; type pinUniverse; radii (0.0); fills (mat5); } + zone6 { id 6; type pinUniverse; radii (0.0); fills (mat6); } + zone7 { id 7; type pinUniverse; radii (0.0); fills (mat7); } + zone8 { id 8; type pinUniverse; radii (0.0); fills (mat8); } + zone9 { id 9; type pinUniverse; radii (0.0); fills (mat9); } + zone10 { id 10; type pinUniverse; radii (0.0); fills (mat10); } + zone11 { id 11; type pinUniverse; radii (0.0); fills (mat11); } + zone12 { id 12; type pinUniverse; radii (0.0); fills (mat12); } + zone13 { id 13; type pinUniverse; radii (0.0); fills (mat13); } + zone14 { id 14; type pinUniverse; radii (0.0); fills (mat14); } + zone15 { id 15; type pinUniverse; radii (0.0); fills (mat15); } + zone16 { id 16; type pinUniverse; radii (0.0); fills (mat16); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat2 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat3 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat4 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat5 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat6 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat7 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat8 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat9 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat10 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat11 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat12 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat13 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat14 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat15 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + mat16 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.25; } + + } + +} + + + diff --git a/InputFiles/IMC/oldInputs/marshakWave32 b/InputFiles/IMC/oldInputs/marshakWave32 new file mode 100644 index 000000000..cfe5c7527 --- /dev/null +++ b/InputFiles/IMC/oldInputs/marshakWave32 @@ -0,0 +1,139 @@ +// Marshak wave simulation using 32 equal spatial regions +// *** lightSpeed and radiationConstant must both be changed to ONE in universalVariables *** + +type IMCPhysicsPackage; + +pop 5000; +limit 50000; +steps 1000; +timeStepSize 0.5; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +source { + type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 1000; +} + +tally { + } + +//grid { dimensions (10 1 1); searchN (100 1 1); } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); } + } + + cells + { + } + universes + { + root { id 100; type rootUniverse; border 1; fill u<200>; } + + lat { id 200; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (0.125 1.0 1.0); + shape (32 1 1); + padMat void; + map ( 1 2 3 4 5 6 7 8 9 10 + 11 12 13 14 15 16 17 18 19 20 + 21 22 23 24 25 26 27 28 29 30 + 31 32); + } + + zone1 { id 1; type pinUniverse; radii (0.0); fills (mat1); } + zone2 { id 2; type pinUniverse; radii (0.0); fills (mat2); } + zone3 { id 3; type pinUniverse; radii (0.0); fills (mat3); } + zone4 { id 4; type pinUniverse; radii (0.0); fills (mat4); } + zone5 { id 5; type pinUniverse; radii (0.0); fills (mat5); } + zone6 { id 6; type pinUniverse; radii (0.0); fills (mat6); } + zone7 { id 7; type pinUniverse; radii (0.0); fills (mat7); } + zone8 { id 8; type pinUniverse; radii (0.0); fills (mat8); } + zone9 { id 9; type pinUniverse; radii (0.0); fills (mat9); } + zone10 { id 10; type pinUniverse; radii (0.0); fills (mat10); } + zone11 { id 11; type pinUniverse; radii (0.0); fills (mat11); } + zone12 { id 12; type pinUniverse; radii (0.0); fills (mat12); } + zone13 { id 13; type pinUniverse; radii (0.0); fills (mat13); } + zone14 { id 14; type pinUniverse; radii (0.0); fills (mat14); } + zone15 { id 15; type pinUniverse; radii (0.0); fills (mat15); } + zone16 { id 16; type pinUniverse; radii (0.0); fills (mat16); } + zone17 { id 17; type pinUniverse; radii (0.0); fills (mat17); } + zone18 { id 18; type pinUniverse; radii (0.0); fills (mat18); } + zone19 { id 19; type pinUniverse; radii (0.0); fills (mat19); } + zone20 { id 20; type pinUniverse; radii (0.0); fills (mat20); } + zone21 { id 21; type pinUniverse; radii (0.0); fills (mat21); } + zone22 { id 22; type pinUniverse; radii (0.0); fills (mat22); } + zone23 { id 23; type pinUniverse; radii (0.0); fills (mat23); } + zone24 { id 24; type pinUniverse; radii (0.0); fills (mat24); } + zone25 { id 25; type pinUniverse; radii (0.0); fills (mat25); } + zone26 { id 26; type pinUniverse; radii (0.0); fills (mat26); } + zone27 { id 27; type pinUniverse; radii (0.0); fills (mat27); } + zone28 { id 28; type pinUniverse; radii (0.0); fills (mat28); } + zone29 { id 29; type pinUniverse; radii (0.0); fills (mat29); } + zone30 { id 30; type pinUniverse; radii (0.0); fills (mat30); } + zone31 { id 31; type pinUniverse; radii (0.0); fills (mat31); } + zone32 { id 32; type pinUniverse; radii (0.0); fills (mat32); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat2 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat3 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat4 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat5 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat6 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat7 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat8 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat9 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat10 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat11 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat12 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat13 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat14 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat15 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat16 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat17 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat18 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat19 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat20 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat21 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat22 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat23 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat24 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat25 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat26 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat27 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat28 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat29 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat30 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat31 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + mat32 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.125; } + + } + +} + + + diff --git a/InputFiles/IMC/oldInputs/marshakWave64 b/InputFiles/IMC/oldInputs/marshakWave64 new file mode 100644 index 000000000..3e8e7ee29 --- /dev/null +++ b/InputFiles/IMC/oldInputs/marshakWave64 @@ -0,0 +1,208 @@ +// Marshak wave simulation using 64 equal spatial regions +// *** lightSpeed and radiationConstant must both be changed to ONE in universalVariables *** + +type IMCPhysicsPackage; + +pop 64000; +limit 2560000; +steps 1000; +timeStepSize 0.5; +printUpdates 4; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; cutoff 0.0; + } + +source { + type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 1000; +} + +tally { + } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); } + } + + cells + { + } + universes + { + root { id 100; type rootUniverse; border 1; fill u<200>; } + + lat { id 200; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (0.0625 1.0 1.0); + shape (64 1 1); + padMat void; + map ( 1 2 3 4 5 6 7 8 9 10 + 11 12 13 14 15 16 17 18 19 20 + 21 22 23 24 25 26 27 28 29 30 + 31 32 33 34 35 36 37 38 39 40 + 41 42 43 44 45 46 47 48 49 50 + 51 52 53 54 55 56 57 58 59 60 + 61 62 63 64); + } + + zone1 { id 1; type pinUniverse; radii (0.0); fills (mat1); } + zone2 { id 2; type pinUniverse; radii (0.0); fills (mat2); } + zone3 { id 3; type pinUniverse; radii (0.0); fills (mat3); } + zone4 { id 4; type pinUniverse; radii (0.0); fills (mat4); } + zone5 { id 5; type pinUniverse; radii (0.0); fills (mat5); } + zone6 { id 6; type pinUniverse; radii (0.0); fills (mat6); } + zone7 { id 7; type pinUniverse; radii (0.0); fills (mat7); } + zone8 { id 8; type pinUniverse; radii (0.0); fills (mat8); } + zone9 { id 9; type pinUniverse; radii (0.0); fills (mat9); } + zone10 { id 10; type pinUniverse; radii (0.0); fills (mat10); } + zone11 { id 11; type pinUniverse; radii (0.0); fills (mat11); } + zone12 { id 12; type pinUniverse; radii (0.0); fills (mat12); } + zone13 { id 13; type pinUniverse; radii (0.0); fills (mat13); } + zone14 { id 14; type pinUniverse; radii (0.0); fills (mat14); } + zone15 { id 15; type pinUniverse; radii (0.0); fills (mat15); } + zone16 { id 16; type pinUniverse; radii (0.0); fills (mat16); } + zone17 { id 17; type pinUniverse; radii (0.0); fills (mat17); } + zone18 { id 18; type pinUniverse; radii (0.0); fills (mat18); } + zone19 { id 19; type pinUniverse; radii (0.0); fills (mat19); } + zone20 { id 20; type pinUniverse; radii (0.0); fills (mat20); } + zone21 { id 21; type pinUniverse; radii (0.0); fills (mat21); } + zone22 { id 22; type pinUniverse; radii (0.0); fills (mat22); } + zone23 { id 23; type pinUniverse; radii (0.0); fills (mat23); } + zone24 { id 24; type pinUniverse; radii (0.0); fills (mat24); } + zone25 { id 25; type pinUniverse; radii (0.0); fills (mat25); } + zone26 { id 26; type pinUniverse; radii (0.0); fills (mat26); } + zone27 { id 27; type pinUniverse; radii (0.0); fills (mat27); } + zone28 { id 28; type pinUniverse; radii (0.0); fills (mat28); } + zone29 { id 29; type pinUniverse; radii (0.0); fills (mat29); } + zone30 { id 30; type pinUniverse; radii (0.0); fills (mat30); } + zone31 { id 31; type pinUniverse; radii (0.0); fills (mat31); } + zone32 { id 32; type pinUniverse; radii (0.0); fills (mat32); } + + zone33 { id 33; type pinUniverse; radii (0.0); fills (mat33); } + zone34 { id 34; type pinUniverse; radii (0.0); fills (mat34); } + zone35 { id 35; type pinUniverse; radii (0.0); fills (mat35); } + zone36 { id 36; type pinUniverse; radii (0.0); fills (mat36); } + zone37 { id 37; type pinUniverse; radii (0.0); fills (mat37); } + zone38 { id 38; type pinUniverse; radii (0.0); fills (mat38); } + zone39 { id 39; type pinUniverse; radii (0.0); fills (mat39); } + zone40 { id 40; type pinUniverse; radii (0.0); fills (mat40); } + zone41 { id 41; type pinUniverse; radii (0.0); fills (mat41); } + zone42 { id 42; type pinUniverse; radii (0.0); fills (mat42); } + zone43 { id 43; type pinUniverse; radii (0.0); fills (mat43); } + zone44 { id 44; type pinUniverse; radii (0.0); fills (mat44); } + zone45 { id 45; type pinUniverse; radii (0.0); fills (mat45); } + zone46 { id 46; type pinUniverse; radii (0.0); fills (mat46); } + zone47 { id 47; type pinUniverse; radii (0.0); fills (mat47); } + zone48 { id 48; type pinUniverse; radii (0.0); fills (mat48); } + zone49 { id 49; type pinUniverse; radii (0.0); fills (mat49); } + zone50 { id 50; type pinUniverse; radii (0.0); fills (mat50); } + zone51 { id 51; type pinUniverse; radii (0.0); fills (mat51); } + zone52 { id 52; type pinUniverse; radii (0.0); fills (mat52); } + zone53 { id 53; type pinUniverse; radii (0.0); fills (mat53); } + zone54 { id 54; type pinUniverse; radii (0.0); fills (mat54); } + zone55 { id 55; type pinUniverse; radii (0.0); fills (mat55); } + zone56 { id 56; type pinUniverse; radii (0.0); fills (mat56); } + zone57 { id 57; type pinUniverse; radii (0.0); fills (mat57); } + zone58 { id 58; type pinUniverse; radii (0.0); fills (mat58); } + zone59 { id 59; type pinUniverse; radii (0.0); fills (mat59); } + zone60 { id 60; type pinUniverse; radii (0.0); fills (mat60); } + zone61 { id 61; type pinUniverse; radii (0.0); fills (mat61); } + zone62 { id 62; type pinUniverse; radii (0.0); fills (mat62); } + zone63 { id 63; type pinUniverse; radii (0.0); fills (mat63); } + zone64 { id 64; type pinUniverse; radii (0.0); fills (mat64); } + + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat2 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat3 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat4 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat5 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat6 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat7 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat8 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat9 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat10 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat11 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat12 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat13 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat14 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat15 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat16 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat17 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat18 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat19 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat20 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat21 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat22 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat23 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat24 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat25 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat26 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat27 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat28 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat29 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat30 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat31 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat32 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + + mat33 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat34 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat35 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat36 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat37 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat38 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat39 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat40 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat41 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat42 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat43 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat44 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat45 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat46 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat47 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat48 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat49 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat50 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat51 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat52 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat53 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat54 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat55 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat56 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat57 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat58 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat59 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat60 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat61 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat62 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat63 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + mat64 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.0625; } + + } + +} + + + diff --git a/InputFiles/IMC/oldInputs/marshakWave8 b/InputFiles/IMC/oldInputs/marshakWave8 new file mode 100644 index 000000000..76aaf57e1 --- /dev/null +++ b/InputFiles/IMC/oldInputs/marshakWave8 @@ -0,0 +1,88 @@ +// Marshak wave simulation using 8 equal spatial regions +// *** lightSpeed and radiationConstant must both be changed to ONE in universalVariables *** + +type IMCPhysicsPackage; + +pop 500; +limit 5000; +steps 10000; +timeStepSize 0.05; +printUpdates 8; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; + } + +source { + type bbSurfaceSource; r (-2 -2 -0.5 0.5 -0.5 0.5); temp 1; N 1000; +} + +tally { + } + +geometry { + type geometryStd; + boundary (0 0 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2 0.5 0.5); } + } + + cells + { + } + universes + { + root { id 100; type rootUniverse; border 1; fill u<200>; } + + lat { id 200; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (0.5 1.0 1.0); + shape (8 1 1); + padMat void; + map ( 1 2 3 4 5 6 7 8); + + } + + zone1 { id 1; type pinUniverse; radii (0.0); fills (mat1); } + zone2 { id 2; type pinUniverse; radii (0.0); fills (mat2); } + zone3 { id 3; type pinUniverse; radii (0.0); fills (mat3); } + zone4 { id 4; type pinUniverse; radii (0.0); fills (mat4); } + zone5 { id 5; type pinUniverse; radii (0.0); fills (mat5); } + zone6 { id 6; type pinUniverse; radii (0.0); fills (mat6); } + zone7 { id 7; type pinUniverse; radii (0.0); fills (mat7); } + zone8 { id 8; type pinUniverse; radii (0.0); fills (mat8); } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat2 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat3 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat4 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat5 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat6 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat7 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + mat8 { temp 0.01; composition {} xsFile ./dataFiles/imcData; volume 0.5; } + + } + +} + + + diff --git a/InputFiles/IMC/sampleInput b/InputFiles/IMC/sampleInput new file mode 100644 index 000000000..7342876a5 --- /dev/null +++ b/InputFiles/IMC/sampleInput @@ -0,0 +1,103 @@ +// +// Intended as a sample/tutorial input file for calculations using IMC physics +// package, to detail settings that differ to other calculation types +// +// Note that a lot of the input files in this branch require lightSpeed and radiationConstant both +// to be changed to ONE in universalVariables before compiling + +type IMCPhysicsPackage; + +pop 1000; + // Maximum total number of particles to be emitted during each time step from all material. + // This number is split between material regions based on the energy they are emitting and is + // reduced if limit is going to be reached. + +limit 10000; + // Sets the maximum size of particle dungeons. Runtime is very dependent on this value so should + // not be set arbitrarily large. + +steps 50; + // The number of time steps to be used in the calculation + +timeStepSize 0.1; + // The time step size for the calculation in seconds + +printUpdates 1; + // The number maximum number of material updates to print to screen. If 0, no updates will be + // printed. + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTimeHT; cutoff 0.7; + } + +// No tallies are required for calculation, but empty dictionary must be given +tally { + } + + +// Geometry is as in all other calculation types. +// Here a simple infinite region is given (a perfectly reflected 1x1x1 cube) + +geometry { + type geometryStd; + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces { + squareBound { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 0.5 0.5 0.5); } + } + + cells { + } + + universes { + root { id 1; type rootUniverse; border 1; fill mat; } + } +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + // Dictionary containing all materials used in geometry + // If desired to have spatial temperature variation, split geometry (above) into desired cells + // and set each cell fill as a DIFFERENT material (e.g. mat1, mat2, mat3, ...) then define + // all materials here. Even if each each mat input is identical, a unique material object + // will be created allowing for a unique temperature evolution. The same xsFile may be used + // for different materials if desired. + + materials { + + // Example: mat + mat { + + temp 1; + // Initial temperature of material [keV]. + + composition {} + // Empty dictionary required for composition. + + xsFile ./DataFiles/sampleData; + // Location of material data file containing material properties. + + volume 1; + // Total volume that this material occupies [cm3], for now need to calculate by hand + // and enter here. May be room to make this automatic in the future. + + } + + // Example 2: mat2 + //mat2 { temp 1; composition {} xsFile ./DataFiles/sampleData2; volume 1 } + +} + +} + + + diff --git a/NuclearData/CMakeLists.txt b/NuclearData/CMakeLists.txt index 47126f0d0..2bdfee55e 100644 --- a/NuclearData/CMakeLists.txt +++ b/NuclearData/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(ceNeutronData) add_subdirectory(mgNeutronData) +add_subdirectory(mgIMCData) add_subdirectory(testNeutronData) add_subdirectory(xsPackages) add_subdirectory(emissionENDF) @@ -11,6 +12,7 @@ add_sources(./nuclearDatabase_inter.f90 ./materialMenu_mod.f90 ./nuclearDataReg_mod.f90 ./neutronMaterial_inter.f90 + ./IMCMaterial_inter.f90 ./Reactions/reactionHandle_inter.f90 ./Reactions/uncorrelatedReactionCE_inter.f90 ./Reactions/reactionMG_inter.f90 diff --git a/NuclearData/IMCMaterial_inter.f90 b/NuclearData/IMCMaterial_inter.f90 new file mode 100644 index 000000000..7864083b4 --- /dev/null +++ b/NuclearData/IMCMaterial_inter.f90 @@ -0,0 +1,146 @@ +module IMCMaterial_inter + + use numPrecision + use particle_class, only : particle + use RNG_class, only : RNG + + ! Nuclear Data Interfaces + use materialHandle_inter, only : materialHandle + use IMCXsPackages_class, only : IMCMacroXSs + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: IMCMaterial_CptrCast + + !! + !! Abstract interface far all IMC Materials (CE and MG) + !! + !! It was created to expose access to some key information in the context of + !! tallying where one is not interested whether MG or CE data is used + !! + !! Interface: + !! materialHandle interface + !! getMacroXSs -> Return Macroscopic XSs given particle with energy data + !! + type, public, abstract, extends(materialHandle) :: IMCMaterial + private + contains + generic :: getMacroXSs => getMacroXSs_byP + procedure(getMacroXSs_byP), deferred :: getMacroXSs_byP + procedure(updateMat), deferred :: updateMat + procedure(getEmittedRad), deferred :: getEmittedRad + procedure(getFleck), deferred :: getFleck + procedure(getTemp), deferred :: getTemp + procedure(setTimeStep), deferred :: setTimeStep + end type IMCMaterial + + abstract interface + + !! + !! Return Macroscopic XSs for the material given particle + !! + !! Args: + !! xss [out] -> Cross section package to store the data + !! p [in] -> Particle that provides energy or energy group + !! + !! Errors: + !! fatalError if energy value/group is outside bounds + !! fatalError if MG particle is given to CE data and vice versa + !! + subroutine getMacroXSs_byP(self, xss, p) + import :: IMCMaterial, particle, IMCMacroXSs + class(IMCMaterial), intent(in) :: self + type(IMCMacroXSs), intent(out) :: xss + class(particle), intent(in) :: p + end subroutine getMacroXSs_byP + + !! + !! Update material properties at each time step + !! First update energy using simple balance, then solve for temperature, + !! then update temperature-dependent properties + !! + !! Args: + !! tallyEnergy [in] -> Energy absorbed into material + !! printUpdate [in, optional] -> Bool, if true then will print updates to screen + !! + subroutine updateMat(self, tallyEnergy, printUpdate) + import :: IMCMaterial, defReal, defBool + class(IMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + logical(defBool), intent(in), optional :: printUpdate + end subroutine updateMat + + !! + !! Return the equilibrium radiation energy density, U_r + !! + function getEmittedRad(self) result(emittedRad) + import :: IMCMaterial, defReal, RNG + class(IMCMaterial), intent(inout) :: self + real(defReal) :: emittedRad + end function getEmittedRad + + !! + !! Get Fleck factor of material + !! + function getFleck(self) result(fleck) + import :: IMCMaterial, defReal + class(IMCMaterial), intent(in) :: self + real(defReal) :: fleck + end function getFleck + + !! + !! Get temperature of material + !! + function getTemp(self) result(T) + import :: IMCMaterial, defReal + class(IMCMaterial), intent(inout) :: self + real(defReal) :: T + end function getTemp + + !! + !! Provide material with time step size + !! + !! Args: + !! dt [in] -> time step size [s] + !! + subroutine setTimeStep(self, dt) + import :: IMCMaterial, defReal + class(IMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: dt + end subroutine setTimeStep + + end interface + +contains + + + !! + !! Cast materialHandle pointer to IMCMaterial pointer + !! + !! Args: + !! source [in] -> source pointer of class materialHandle + !! + !! Result: + !! Null if source is not of IMCMaterial + !! Pointer to source if source is IMCMaterial class + !! + pure function IMCMaterial_CptrCast(source) result(ptr) + class(materialHandle), pointer, intent(in) :: source + class(IMCMaterial), pointer :: ptr + + select type(source) + class is(IMCMaterial) + ptr => source + + class default + ptr => null() + end select + + end function IMCMaterial_CptrCast + + +end module IMCMaterial_inter diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 3236083d2..36d69f7b0 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -10,13 +10,16 @@ !! nameMap -> Map that maps material name to matIdx !! !! Interface: -!! init -> Load material definitions from a dictionary -!! kill -> Return to uninitialised state -!! display -> Display information about all defined materials to console -!! getMatPtr -> Return pointer to a detailed material information (materialItem) -!! nMat -> Return number of materials -!! matName -> Return material Name given Index -!! matIdx -> Return material Index given Name +!! init -> Load material definitions from a dictionary +!! kill -> Return to uninitialised state +!! display -> Display information about all defined materials to console +!! getMatPtr -> Return pointer to a detailed material information (materialItem) +!! nMat -> Return number of materials +!! matName -> Return material Name given Index +!! matTemp -> Return material Temperature given Index +!! matVol -> Return material Volume given Index +!! matFile -> Return file path to material data given matIdx +!! matIdx -> Return material Index given Name !! module materialMenu_mod @@ -64,6 +67,7 @@ module materialMenu_mod !! name -> name of material !! matIdx -> material index of the material !! T -> material temperature [K] + !! V -> volume of material zone [cm3] !! dens -> vector of densities [1/barn/cm] !! nuclides -> associated vector of nuclide types !! extraInfo -> dictionary with extra keywords @@ -93,6 +97,7 @@ module materialMenu_mod character(nameLen) :: name = '' integer(shortInt) :: matIdx = 0 real(defReal) :: T = ZERO + real(defReal) :: V = ZERO real(defReal),dimension(:),allocatable :: dens type(nuclideInfo),dimension(:),allocatable :: nuclides type(dictionary) :: extraInfo @@ -112,6 +117,9 @@ module materialMenu_mod public :: getMatPtr public :: nMat public :: matName + public :: matTemp + public :: matVol + public :: matFile public :: matIdx contains @@ -131,7 +139,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 @@ -204,7 +212,7 @@ end subroutine display !! Result: !! nameLen long character with material name !! - !! Erorrs: + !! Errors: !! If idx is -ve or larger then number of defined materials !! Empty string '' is returned as its name !! @@ -221,6 +229,73 @@ function matName(idx) result(name) end function matName + !! + !! Return starting temperature of material given index + !! + !! Args: + !! idx [in] -> Material Index + !! + !! Result: + !! Temperature of material as given in input file + !! + !! Errors: + !! If idx is -ve or larger then number of defined materials + !! then -1 is returned as its temperature + !! + function matTemp(idx) result(T) + integer(shortInt), intent(in) :: idx + real(defReal) :: T + + if(idx <= 0 .or. nMat() < idx) then + T = -ONE + else + T = materialDefs(idx) % T + end if + + end function matTemp + + !! + !! Return volume of material given index + !! + !! Args: + !! idx [in] -> Material Index + !! + !! Result: + !! Volume of material as given in input file + !! + !! Errors: + !! If idx is -ve or larger then number of defined materials + !! then -1 is returned as its volume + !! + function matVol(idx) result(vol) + integer(shortInt), intent(in) :: idx + real(defReal) :: vol + + if( idx <= 0 .or. nMat() < idx) then + vol = -ONE + else + vol = materialDefs(idx) % V + end if + + end function matVol + + !! + !! Return file path to material XS data given index + !! + !! Args: + !! idx [in] -> Material index + !! + !! Result: + !! Path to file + !! + function matFile(idx) result(path) + integer(shortInt), intent(in) :: idx + character(pathLen) :: path + + call materialDefs(idx) % extraInfo % get(path,'xsFile') + + end function matFile + !! !! Return material index Given Name !! @@ -270,6 +345,7 @@ subroutine init_materialItem(self, name, dict) ! Load easy components c self % name = name call dict % get(self % T,'temp') + call dict % getOrDefault(self % V, 'volume', ZERO) ! Get composition dictionary and load composition compDict => dict % getDictPtr('composition') @@ -318,6 +394,7 @@ subroutine kill_materialItem(self) self % name = '' self % matIdx = 0 self % T = ZERO + self % V = ZERO ! Deallocate allocatable components if(allocated(self % dens)) deallocate(self % dens) diff --git a/NuclearData/mgIMCData/CMakeLists.txt b/NuclearData/mgIMCData/CMakeLists.txt new file mode 100644 index 000000000..dd219aafc --- /dev/null +++ b/NuclearData/mgIMCData/CMakeLists.txt @@ -0,0 +1,8 @@ +# Add source files for compilation +add_sources(./mgIMCMaterial_inter.f90 + ./mgIMCDatabase_inter.f90 + ./baseMgIMC/baseMgIMCMaterial_class.f90 + ./baseMgIMC/baseMgIMCDatabase_class.f90) + +# Add tests +#add_integration_tests(./baseMgIMC/Tests/baseMgIMCDatabase_iTest.f90) diff --git a/NuclearData/mgIMCData/baseMgIMC/baseMgIMCDatabase_class.f90 b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCDatabase_class.f90 new file mode 100644 index 000000000..30cc86171 --- /dev/null +++ b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCDatabase_class.f90 @@ -0,0 +1,438 @@ +module baseMgIMCDatabase_class + + use numPrecision + use endfConstants + use universalVariables, only : VOID_MAT + use genericProcedures, only : fatalError, numToChar + use particle_class, only : particle + use charMap_class, only : charMap + use dictionary_class, only : dictionary + use dictParser_func, only : fileToDict + + ! Nuclear Data Interfaces + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase + use materialHandle_inter, only : materialHandle + use nuclideHandle_inter, only : nuclideHandle + use reactionHandle_inter, only : reactionHandle + use materialMenu_mod, only : materialItem, mm_getMatPtr => getMatPtr, mm_nMat => nMat, & + mm_nameMap => nameMap + + ! baseMgIMC Objects + use baseMgIMCMaterial_class, only : baseMgIMCMaterial + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: baseMgIMCDatabase_TptrCast + public :: baseMgIMCDatabase_CptrCast + + !! + !! Basic type of MG nuclear Data for IMCs + !! + !! All materials in aproblem are baseMgMaterials. See its documentation for + !! details on how the physics is handled + !! + !! Public Members: + !! mats -> array containing all defined materials (by matIdx) + !! activeMats -> list of matIdxs of materials active in the problem + !! + !! Interface: + !! nuclearDatabase interface + !! + type, public, extends(mgIMCDatabase) :: baseMgIMCDatabase + type(baseMgIMCMaterial), dimension(:), pointer :: mats => null() + integer(shortInt), dimension(:), allocatable :: activeMats + integer(shortInt) :: nG = 0 + + contains + ! Superclass Interface + procedure :: getTransMatXS + procedure :: getTotalMatXS + procedure :: getMajorantXS + procedure :: matNamesMap + procedure :: getMaterial + procedure :: getNuclide + procedure :: getReaction + procedure :: getEmittedRad + procedure :: updateProperties + procedure :: setTimeStep + procedure :: kill + procedure :: init + procedure :: activate + + ! Local interface + procedure :: nGroups + + end type baseMgIMCDatabase + +contains + + !! + !! Get Transport XS given a particle + !! + !! See nuclearDatabase documentation for details + !! + !! Note: + !! DOES NOT check if particle is MG. Will refer to G in the particle and give error + !! if the value is invalid + !! + function getTransMatXS(self, p, matIdx) result(xs) + class(baseMgIMCDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + + xs = self % getTotalMatXS(p, matIdx) + + end function getTransMatXS + + !! + !! Get Total XS given a particle + !! + !! See nuclearDatabase documentation for details + !! + !! Note: + !! DOES NOT check if particle is MG. Will refer to G in the particle and give error + !! if the value is invalid + !! + function getTotalMatXS(self, p, matIdx) result(xs) + class(baseMgIMCDatabase), intent(inout) :: self + class(particle), intent(in) :: p + integer(shortInt), intent(in) :: matIdx + real(defReal) :: xs + + if (matIdx == VOID_MAT) then + xs = ZERO + else + xs = self % mats(matIdx) % getTotalXS(p % G, p % pRNG) + end if + + end function getTotalMatXS + + !! + !! Get Majorant XS given a particle + !! + !! See nuclearDatabase documentation for details + !! + !! Note: + !! DOES NOT check if particle is MG. Will refer to G in the particle and give error + !! if the value is invalid + !! + function getMajorantXS(self, p) result(xs) + class(baseMgIMCDatabase), intent(inout) :: self + class(particle), intent(in) :: p + real(defReal) :: xs + integer(shortInt) :: i, idx + + xs = ZERO + do i=1,size(self % activeMats) + idx = self % activeMats(i) + xs = max(xs, self % getTotalMatXS(p, idx)) + end do + + end function getMajorantXS + + !! + !! Return pointer to material names map + !! + !! See nuclearDatabase documentation for details + !! + function matNamesMap(self) result(map) + class(baseMgIMCDatabase), intent(in) :: self + type(charMap), pointer :: map + + map => mm_nameMap + + end function matNamesMap + + !! + !! Return pointer to a material in the database + !! + !! See nuclearDatabase documentation for details + !! + function getMaterial(self, matIdx) result(mat) + class(baseMgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in) :: matIdx + class(materialHandle), pointer :: mat + + if(matIdx < 1 .or. matIdx > size(self % mats)) then + mat => null() + else + mat => self % mats(matIdx) + end if + + end function getMaterial + + !! + !! Return pointer to a nuclide in the database + !! + !! See nuclearDatabase documentation for details + !! + !! Note: + !! This database has no nucldie. Returns NULL always! + !! + function getNuclide(self, nucIdx) result(nuc) + class(baseMgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in) :: nucIdx + class(nuclideHandle), pointer :: nuc + + nuc => null() + + end function getNuclide + + !! + !! Return pointer to a reaction + !! + !! See nuclearDatabase documentation for details + !! + function getReaction(self, MT, idx) result(reac) + class(baseMgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in) :: MT + integer(shortInt), intent(in) :: idx + class(reactionHandle), pointer :: reac + character(100), parameter :: Here = 'getReaction (baseMgIMCDatabase_class.f90)' + + reac => null() + call fatalError(Here, "Pointless function call") + + end function getReaction + + !! + !! Return energy to be emitted during current time step + !! + !! Args: + !! matIdx [in] [optional] -> If provided, return the energy to be emitted from only matIdx + !! Otherwise, return total energy to be emitted from all mats + !! + function getEmittedRad(self, matIdx) result(energy) + class(baseMgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in), optional :: matIdx + real(defReal) :: energy + integer(shortInt) :: i + + ! If matIdx provided, return radiation emitted from only that material + if (present(matIdx)) then + energy = self % mats(matIdx) % getEmittedRad() + return + end if + + ! Otherwise, return total energy emitted from all materials + energy = 0 + + do i=1, size(self % mats) + energy = energy + self % mats(i) % getEmittedRad() + end do + + end function getEmittedRad + + !! + !! Update material properties based on energy absorbed during the time step + !! + subroutine updateProperties(self, tallyEnergy, printUpdates) + class(baseMgIMCDatabase), intent(inout) :: self + real(defReal), dimension(:), intent(in) :: tallyEnergy + integer(shortInt), intent(in) :: printUpdates + integer(shortInt) :: i + character(100), parameter :: Here = 'updateProperties (baseMgIMCDatabase_class.f90)' + + ! Check for valid inputs + if (size(tallyEnergy) /= size(self % mats)) call fatalError(Here, & + &'Energy tally array must have size nMats') + if (printUpdates > size(self % mats)) call fatalError(Here, & + &'printUpdates must be <= nMats') + + ! Update mats to be printed (if any) + do i = 1, printUpdates + call self % mats(i) % updateMat(tallyEnergy(i), .true.) + end do + + ! Update remaining mats + !$omp parallel do + do i = (printUpdates+1), size(tallyEnergy) + call self % mats(i) % updateMat(tallyEnergy(i)) + end do + !$omp end parallel do + + end subroutine updateProperties + + !! + !! Provide each material with time step to calculate initial fleck factor + !! + subroutine setTimeStep(self, deltaT) + class(baseMgIMCDatabase), intent(inout) :: self + real(defReal), intent(in) :: deltaT + integer(shortInt) :: i + + do i=1, size(self % mats) + call self % mats(i) % setTimeStep(deltaT) + end do + + end subroutine setTimeStep + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(baseMgIMCDatabase), intent(inout) :: self + + if(associated(self % mats)) then + call self % mats % kill() + deallocate(self % mats) + end if + + if(allocated(self % activeMats)) deallocate (self % activeMats) + self % nG = 0 + + end subroutine kill + + !! + !! Initialise Database from dictionary and pointer to self + !! + !! See nuclearDatabase documentation for details + !! + subroutine init(self, dict, ptr, silent) + class(baseMgIMCDatabase), target,intent(inout) :: self + class(dictionary), intent(in) :: dict + class(nuclearDatabase), pointer,intent(in) :: ptr + logical(defBool), intent(in), optional :: silent + logical(defBool) :: loud + integer(shortInt) :: i, nMat + type(materialItem), pointer :: matDef + character(pathLen) :: path + type(dictionary) :: tempDict + character(100), parameter :: Here = 'init (baseMgIMCDatabase_class.f90)' + + ! Prevent reallocations + call self % kill() + + ! Set build console output flag + if(present(silent)) then + loud = .not.silent + else + loud = .true. + end if + + ! Find number of materials and allocate space + nMat = mm_nMat() + + allocate(self % mats(nMat)) + + ! Build materials + do i=1,nMat + + ! Get Path to the xsFile + matDef => mm_getMatPtr(i) + call matDef % extraInfo % get(path,'xsFile') + + ! Print status + if(loud) then + print '(A)', "Building material: " // trim(matDef % name) // " From: " // trim(path) + end if + + ! Load dictionary + call fileToDict(tempDict, path) + + ! Add temperature and volume into dictionary + call tempDict % store('T', matDef % T) + call tempDict % store('V', matdef % V) + + ! Initialise material + call self % mats(i) % init(tempDict) + + end do + + ! Load and verify number of groups + self % nG = self % mats(1) % nGroups() + do i=2,nMat + if(self % nG /= self % mats(i) % nGroups()) then + call fatalError(Here,'Inconsistant # of groups in materials in matIdx'//numToChar(i)) + end if + end do + + end subroutine init + + !! + !! Activate this nuclearDatabase + !! + !! See nuclearDatabase documentation for details + !! + subroutine activate(self, activeMat) + class(baseMgIMCDatabase), intent(inout) :: self + integer(shortInt), dimension(:), intent(in) :: activeMat + + if(allocated(self % activeMats)) deallocate(self % activeMats) + self % activeMats = activeMat + + end subroutine activate + + !! + !! Return number of energy groups in this database + !! + !! Args: + !! None + !! + !! Errors: + !! None + !! + pure function nGroups(self) result(nG) + class(baseMgIMCDatabase), intent(in) :: self + integer(shortInt) :: nG + + nG = self % nG + + end function nGroups + + !! + !! Cast nuclearDatabase pointer to baseMgIMCDatabase type pointer + !! + !! Args: + !! source [in] -> source pointer of class nuclearDatabase + !! + !! Result: + !! Null if source is not of baseMgIMCDatabase type + !! Target points to source if source is baseMgIMCDatabasetype + !! + pure function baseMgIMCDatabase_TptrCast(source) result(ptr) + class(nuclearDatabase), pointer, intent(in) :: source + type(baseMgIMCDatabase), pointer :: ptr + + select type(source) + type is(baseMgIMCDatabase) + ptr => source + + class default + ptr => null() + end select + + end function baseMgIMCDatabase_TptrCast + + !! + !! Cast nuclearDatabase pointer to baseMgIMCDatabase class pointer + !! + !! Args: + !! source [in] -> source pointer of class nuclearDatabase + !! + !! Result: + !! Null if source is not of baseMgIMCDatabase class + !! Target points to source if source is baseMgIMCDatabase class + !! + pure function baseMgIMCDatabase_CptrCast(source) result(ptr) + class(nuclearDatabase), pointer, intent(in) :: source + class(baseMgIMCDatabase), pointer :: ptr + + select type(source) + class is(baseMgIMCDatabase) + ptr => source + + class default + ptr => null() + end select + + end function baseMgIMCDatabase_CptrCast + + +end module baseMgIMCDatabase_class diff --git a/NuclearData/mgIMCData/baseMgIMC/baseMgIMCMaterial_class.f90 b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCMaterial_class.f90 new file mode 100644 index 000000000..3560bc1ff --- /dev/null +++ b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCMaterial_class.f90 @@ -0,0 +1,503 @@ +module baseMgIMCMaterial_class + + use numPrecision + use endfConstants + use universalVariables + use genericProcedures, only : fatalError, numToChar + use RNG_class, only : RNG + use dictionary_class, only : dictionary + use poly_func + + ! Nuclear Data Interfaces + use materialHandle_inter, only : materialHandle + use mgIMCMaterial_inter, only : mgIMCMaterial, kill_super => kill + use IMCXSPackages_class, only : IMCMacroXSs + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: baseMgIMCMaterial_TptrCast + public :: baseMgIMCMaterial_CptrCast + + ! Public data location parameters + ! Use them if accessing data entries directly + integer(shortInt), parameter, public :: TOTAL_XS = 1 + integer(shortInt), parameter, public :: IESCATTER_XS = 2 + integer(shortInt), parameter, public :: CAPTURE_XS = 3 + integer(shortInt), parameter, public :: PLANCK_XS = 4 + + ! IMC Calculation Type + integer(shortInt), parameter, public :: IMC = 1 + integer(shortInt), parameter, public :: ISMC = 2 + + !! + !! Basic type of MG material data + !! + !! Stores MG data in a table. + !! All other scattering reactions are lumped into single multiplicative scattering, + !! which is stored as INELASTIC scatering in macroXSs package! After all it is inelastic in + !! the sense that outgoing group can change. Diffrent types of multiplicative scattering can be + !! build. See doc of "init" procedure for details. + !! + !! Public members: + !! data -> Rank 2 array with all XSs data + !! + !! Interface: + !! materialHandle interface + !! mgIMCMaterial interface + !! init -> initialise Basic MG Material from dictionary and config keyword + !! nGroups -> returns number of energy groups + !! updateMat -> update material properties as required for IMC calculation + !! getEmittedRad -> returns the radiation to be emitted in current timestep + !! getFleck -> returns current material Fleck factor + !! getTemp -> returns current material temperature + !! + !! Note: + !! Order of "data" array is: data(XS_type, Group #) + !! Dictionary with data must contain following entries: + !! -> numberOfGroups + !! + type, public, extends(mgIMCMaterial) :: baseMgIMCMaterial + real(defReal),dimension(:,:), allocatable :: data + real(defReal),dimension(:), allocatable :: cv + real(defReal),dimension(:), allocatable :: updateEqn + real(defReal),dimension(:), allocatable :: absEqn + real(defReal),dimension(:), allocatable :: scattEqn + real(defReal),dimension(:), allocatable :: planckEqn + real(defReal) :: T + real(defReal) :: V + real(defReal) :: fleck + real(defReal) :: alpha + real(defReal) :: deltaT + real(defReal) :: sigmaP + real(defReal) :: matEnergy + integer(shortInt) :: calcType + + contains + ! Superclass procedures + procedure :: kill + procedure :: getMacroXSs_byG + procedure :: getTotalXS + + ! Local procedures + procedure :: init + procedure :: nGroups + procedure :: updateMat + procedure :: getEmittedRad + procedure :: getFleck + procedure :: getTemp + procedure :: setTimeStep + + procedure, private :: updateMatIMC + procedure, private :: updateMatISMC + procedure, private :: tempFromEnergy + procedure, private :: sigmaFromTemp + + end type baseMgIMCMaterial + +contains + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(baseMgIMCMaterial), intent(inout) :: self + + ! Call superclass procedure + call kill_super(self) + + ! Kill local content + if(allocated(self % data)) deallocate(self % data) + + end subroutine kill + + !! + !! Load Macroscopic XSs into the provided package for a given group index G + !! + !! See mgIMCMaterial documentation for more details + !! + subroutine getMacroXSs_byG(self, xss, G, rand) + class(baseMgIMCMaterial), intent(in) :: self + type(IMCMacroXSs), intent(out) :: xss + integer(shortInt), intent(in) :: G + class(RNG), intent(inout) :: rand + character(100), parameter :: Here = ' getMacroXSs (baseMgIMCMaterial_class.f90)' + + ! Verify bounds + if(G < 1 .or. self % nGroups() < G) then + call fatalError(Here,'Invalid group number: '//numToChar(G)// & + ' Data has only: ' // numToChar(self % nGroups())) + end if + + ! Get XSs + xss % total = self % data(TOTAL_XS, G) + xss % elasticScatter = ZERO + xss % inelasticScatter = self % data(IESCATTER_XS, G) + xss % capture = self % data(CAPTURE_XS, G) + xss % planck = self % data(PLANCK_XS, G) + + end subroutine getMacroXSs_byG + + !! + !! Return Total XSs for energy group G + !! + !! See mgIMCMaterial documentationfor details + !! + function getTotalXS(self, G, rand) result(xs) + class(baseMgIMCMaterial), intent(in) :: self + integer(shortInt), intent(in) :: G + class(RNG), intent(inout) :: rand + real(defReal) :: xs + character(100), parameter :: Here = ' getTotalXS (baseMgIMCMaterial_class.f90)' + + ! Verify bounds + if(G < 1 .or. self % nGroups() < G) then + call fatalError(Here,'Invalid group number: '//numToChar(G)// & + ' Data has only: ' // numToChar(self % nGroups())) + xs = ZERO ! Avoid warning + end if + xs = self % data(TOTAL_XS, G) + + end function getTotalXS + + + !! + !! Initialise Base MG IMC Material fromdictionary + !! + !! Args: + !! dict [in] -> Input dictionary with all required XSs + !! + !! Errors: + !! FatalError if data in dictionary is invalid (inconsistant # of groups; + !! -ve entries in P0 XSs) + !! + subroutine init(self, dict) + class(baseMgIMCMaterial), intent(inout) :: self + class(dictionary),target, intent(in) :: dict + integer(shortInt) :: nG, N + real(defReal), dimension(:), allocatable :: temp + character(100), parameter :: Here = 'init (baseMgIMCMaterial_class.f90)' + + ! Read number of groups + call dict % get(nG, 'numberOfGroups') + if(nG < 1) call fatalError(Here,'Number of groups is invalid' // numToChar(nG)) + + ! Allocate space for data + N = 4 + allocate(self % data(N, nG)) + + ! Store alpha setting + call dict % getOrDefault(self % alpha, 'alpha', ONE) + + ! Read opacity equations + call dict % get(temp, 'capture') + self % absEqn = temp + call dict % get(temp, 'scatter') + self % scattEqn = temp + + ! Build planck opacity equation + ! For grey case, sigmaP = sigmaA. Will become more complicated for frequency-dependent case + self % planckEqn = self % absEqn + + ! Read heat capacity equation + call dict % get(temp, 'cv') + self % cv = temp + + ! Build update equation + call poly_integrate(temp) + self % updateEqn = temp + + ! Read initial temperature and volume + call dict % get(self % T, 'T') + call dict % get(self % V, 'V') + + ! Calculate initial opacities and energy + call self % sigmaFromTemp() + self % matEnergy = poly_eval(self % updateEqn, self % T) * self % V + + ! Set calculation type (will support ISMC in the future) + self % calcType = IMC + + end subroutine init + + + !! + !! Provide material with time step size + !! + !! Args: + !! dt [in] -> time step size [s] + !! + !! Errors: + !! fatalError if calculation type is invalid (valid options are IMC or ISMC) + !! + subroutine setTimeStep(self, dt) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: dt + real(defReal) :: beta, eta, zeta + character(100), parameter :: Here = 'setTimeStep (baseMgIMCMaterial_class.f90)' + + self % deltaT = dt + + beta = 4 * radiationConstant * self % T**3 / poly_eval(self % cv, self % T) + + ! Use time step size to calculate fleck factor + if(self % calcType == IMC) then + self % fleck = 1/(1+1*self % sigmaP*lightSpeed*beta*self % deltaT*self % alpha) + + else if(self % calcType == ISMC) then + eta = radiationConstant * self % T**4 / self % matEnergy + zeta = beta - eta + self % fleck = 1/(1 + zeta*self % sigmaP*lightSpeed*self % deltaT) + + else + call fatalError(Here, 'Calculation type invalid or not set') + end if + + + end subroutine setTimeStep + + !! + !! Return number of energy groups + !! + !! Args: + !! None + !! + !! Errors: + !! None + !! + pure function nGroups(self) result(nG) + class(baseMgIMCMaterial), intent(in) :: self + integer(shortInt) :: nG + + if(allocated(self % data)) then + nG = size(self % data,2) + else + nG = 0 + end if + + end function nGroups + + !! + !! Cast materialHandle pointer to baseMgIMCMaterial type pointer + !! + !! Args: + !! source [in] -> source pointer of class materialHandle + !! + !! Result: + !! Null if source is not of baseMgIMCMaterial type + !! Target points to source if source is baseMgIMCMaterialtype + !! + pure function baseMgIMCMaterial_TptrCast(source) result(ptr) + class(materialHandle), pointer, intent(in) :: source + type(baseMgIMCMaterial), pointer :: ptr + + select type(source) + type is(baseMgIMCMaterial) + ptr => source + + class default + ptr => null() + end select + + end function baseMgIMCMaterial_TptrCast + + !! + !! Cast materialHandle pointer to baseMgIMCMaterial class pointer + !! + !! Args: + !! source [in] -> source pointer of class materialHandle + !! + !! Result: + !! Null if source is not of baseMgIMCMaterial class + !! Target points to source if source is baseMgIMCMaterial class + !! + pure function baseMgIMCMaterial_CptrCast(source) result(ptr) + class(materialHandle), pointer, intent(in) :: source + class(baseMgIMCMaterial), pointer :: ptr + + select type(source) + class is(baseMgIMCMaterial) + ptr => source + + class default + ptr => null() + end select + + end function baseMgIMCMaterial_CptrCast + + !! + !! Update material properties at each time step + !! First update energy using simple balance, then solve for temperature, + !! then update temperature-dependent properties + !! + !! Args: + !! tallyEnergy [in] -> Energy absorbed into material + !! printUpdate [in, optional] -> Bool, if true then will print updates to screen + !! + subroutine updateMat(self, tallyEnergy, printUpdate) + class(baseMgIMCMaterial),intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + logical(defBool), intent(in), optional :: printUpdate + character(100), parameter :: Here = "updateMat (baseMgIMCMaterial_class.f90)" + + ! Print current properties + if (present(printUpdate)) then + if (printUpdate .eqv. .True.) then + print *, " T_old = ", self % T + print *, " matEnergy at start of timestep =", self % matEnergy + print *, " emittedRad = ", self % getEmittedRad() + print *, " tallyEnergy = ", tallyEnergy + end if + end if + + select case (self % calcType) + + case(IMC) + call self % updateMatIMC(tallyEnergy) + + case(ISMC) + call self % updateMatISMC(tallyEnergy) + + case default + call fatalError(Here, "Invalid calculation type") + + end select + + ! Print updated properties + if (present(printUpdate)) then + if(printUpdate .eqv. .True.) then + print *, " matEnergy at end of timestep = ", self % matEnergy + print *, " T_new = ", self % T + end if + end if + + end subroutine updateMat + + !! + !! Material update for IMC calculation + !! + subroutine updateMatIMC(self, tallyEnergy) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + real(defReal) :: beta + character(100), parameter :: Here = "updateMatIMC (baseMgIMCMaterial_class.f90)" + + ! Return if no energy change + if (self % getEmittedRad() == tallyEnergy) return + + ! Update material internal energy + self % matEnergy = self % matEnergy - self % getEmittedRad() + tallyEnergy + + ! Update material temperature + self % T = self % tempFromEnergy() + + ! Update sigma + call self % sigmaFromTemp() + + if( self % T < 0 ) then + call fatalError(Here, "Temperature is negative") + end if + + beta = 4 * radiationConstant * self % T**3 / poly_eval(self % cv, self % T) + + self % fleck = 1/(1+1*self % sigmaP*lightSpeed*beta*self % deltaT*self % alpha) + + end subroutine updateMatIMC + + !! + !! Material update for ISMC calculation + !! + subroutine updateMatISMC(self, tallyEnergy) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + real(defReal) :: beta, eta, zeta + + ! Update material internal energy + self % matEnergy = tallyEnergy + + ! Update material temperature + self % T = self % tempFromEnergy() + + ! Update ISMC equivalent of fleck factor + beta = 4*radiationConstant * self % T**3 / poly_eval(self % cv, self % T) + eta = radiationConstant * self % T**4 / self % matEnergy + zeta = beta - eta + self % fleck = 1 / (1 + zeta*self % sigmaP*lightSpeed*self % deltaT) + + ! Update sigma + call self % sigmaFromTemp() + + end subroutine updateMatISMC + + !! + !! Calculate the temperature of material from internal energy + !! + function tempFromEnergy(self) result(T) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: T, energyDens + + energyDens = self % matEnergy / self % V + + if (energyDens == 0) then + T = 0 + else + T = poly_solve(self % updateEqn, self % cv, self % T, energyDens) + end if + + end function tempFromEnergy + + !! + !! Calculate sigma from current temp + !! + subroutine sigmaFromTemp(self) + class(baseMgIMCMaterial), intent(inout) :: self + + self % sigmaP = poly_eval(self % planckEqn, self % T) + + self % data(CAPTURE_XS,:) = poly_eval(self % absEqn, self % T) + self % data(IESCATTER_XS,:) = poly_eval(self % scattEqn, self % T) + self % data(TOTAL_XS,:) = self % data(CAPTURE_XS,:) + self % data(IESCATTER_XS,:) + self % data(PLANCK_XS,:) = poly_eval(self % planckEqn, self % T) + + end subroutine sigmaFromTemp + + + !! + !! Return the energy to be emitted during time step, E_r + !! + function getEmittedRad(self) result(emittedRad) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: U_r, emittedRad + + U_r = radiationConstant * (self % T)**4 + + emittedRad = lightSpeed * self % deltaT * self % sigmaP * self % fleck * U_r * self % V + + end function getEmittedRad + + !! + !! Return the fleck factor of the material + !! + function getFleck(self) result(fleck) + class(baseMgIMCMaterial),intent(in) :: self + real(defReal) :: fleck + + fleck = self % fleck + + end function getFleck + + !! + !! Get temperature of material + !! + function getTemp(self) result(T) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: T + + T = self % T + + end function getTemp + +end module baseMgIMCMaterial_class diff --git a/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 b/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 new file mode 100644 index 000000000..e771c51f6 --- /dev/null +++ b/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 @@ -0,0 +1,95 @@ +module mgIMCDatabase_inter + + use numPrecision + + ! Nuclear Data Interfaces & Objects + use nuclearDatabase_inter, only : nuclearDatabase + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: mgIMCDatabase_CptrCast + + !! + !! An abstract class that groups all MG IMC Data objects + !! + !! It does nothing, It adds nothing, + !! It just provides a common superclass for related classes + !! + type, public, abstract, extends(nuclearDatabase) :: mgIMCDatabase + + contains + procedure(getEmittedRad), deferred :: getEmittedRad + procedure(updateProperties), deferred :: updateProperties + procedure(setTimeStep), deferred :: setTimeStep + + end type mgIMCDatabase + + abstract interface + + !! + !! Return energy to be emitted during current time step + !! + !! Args: + !! matIdx [in] [optional] -> If provided, return the energy to be emitted from only matIdx + !! Otherwise, return total energy to be emitted from all mats + !! + function getEmittedRad(self, matIdx) result(energy) + import :: mgIMCDatabase, shortInt, defReal + class(mgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in), optional :: matIdx + real(defReal) :: energy + end function getEmittedRad + + !! + !! Update material properties based on energy absorbed during the time step + !! + subroutine updateProperties(self, tallyEnergy, printUpdates) + import mgIMCDatabase, defReal, shortInt + class(mgIMCDatabase), intent(inout) :: self + real(defReal), dimension(:), intent(in) :: tallyEnergy + integer(shortInt), intent(in) :: printUpdates + end subroutine updateProperties + + !! + !! Provide each material with time step to calculate initial fleck factor + !! + subroutine setTimeStep(self, deltaT) + import mgIMCDatabase, defReal + class(mgIMCDatabase), intent(inout) :: self + real(defReal), intent(in) :: deltaT + end subroutine setTimeStep + + end interface + +contains + + !! + !! Cast nuclearDatabase pointer to mgIMCDatabase class pointer + !! + !! Args: + !! source [in] -> source pointer of class nuclearDatabase + !! + !! Result: + !! Null if source is not of mgIMCDatabase class + !! Target points to source if source is mgIMCDatabase class + !! + pure function mgIMCDatabase_CptrCast(source) result(ptr) + class(nuclearDatabase), pointer, intent(in) :: source + class(mgIMCDatabase), pointer :: ptr + + select type(source) + class is(mgIMCDatabase) + ptr => source + + class default + ptr => null() + end select + + end function mgIMCDatabase_CptrCast + + +end module mgIMCDatabase_inter diff --git a/NuclearData/mgIMCData/mgIMCMaterial_inter.f90 b/NuclearData/mgIMCData/mgIMCMaterial_inter.f90 new file mode 100644 index 000000000..f167cf805 --- /dev/null +++ b/NuclearData/mgIMCData/mgIMCMaterial_inter.f90 @@ -0,0 +1,210 @@ +module mgIMCMaterial_inter + + use numPrecision + use genericProcedures, only : fatalError + use RNG_class, only : RNG + use particle_class, only : particle + + ! Nuclear Data Handles + use materialHandle_inter, only : materialHandle + use IMCMaterial_inter, only : IMCMaterial + use IMCXsPackages_class, only : IMCMacroXSs + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: mgIMCMaterial_CptrCast + + !! + !! Extendable procedures is subclasses + !! + public :: kill + + !! + !! Abstract interface for all MG IMC Materials + !! + !! Interface: + !! materialHandle interface + !! neutroNMaterial interface + !! getMacroXSs -> Get macroscopic XSs directly from group number and RNG + !! + type, public, abstract, extends(IMCMaterial) :: mgIMCMaterial + private + + contains + ! Superclass procedures + procedure :: kill + generic :: getMacroXSs => getMacroXSs_byG + procedure :: getMacroXSs_byP + + ! Local procedures + procedure(getMacroXSs_byG), deferred :: getMacroXSs_byG + procedure(getTotalXS), deferred :: getTotalXS + procedure(updateMat), deferred :: updateMat + procedure(getEmittedRad), deferred :: getEmittedRad + procedure(getFleck), deferred :: getFleck + procedure(getTemp), deferred :: getTemp + procedure(setTimeStep), deferred :: setTimeStep + + end type mgIMCMaterial + + + + abstract interface + + !! + !! Return Macroscopic XSs for the material + !! + !! Args: + !! xss [out] -> Cross section package to store the data + !! G [in] -> Requested energy group + !! rand [inout] -> Random Number Generator + !! + !! Errors: + !! fatalError if G is out-of-bounds for the stored data + !! + subroutine getMacroXSs_byG(self, xss, G, rand) + import :: mgIMCMaterial, IMCMacroXSs, shortInt, RNG + class(mgIMCMaterial), intent(in) :: self + type(IMCMacroXSs), intent(out) :: xss + integer(shortInt), intent(in) :: G + class(RNG), intent(inout) :: rand + end subroutine getMacroXSs_byG + + !! + !! Return Macroscopic Total XS for the material + !! + !! Args: + !! G [in] -> Requested energygroup + !! rand [inout] -> Random number generator + !! + !! Errors: + !! fatalError if G is out-of-bounds for the stored data + !! + function getTotalXS(self, G, rand) result(xs) + import :: mgIMCMaterial, defReal, shortInt, RNG + class(mgIMCMaterial), intent(in) :: self + integer(shortInt), intent(in) :: G + class(RNG), intent(inout) :: rand + real(defReal) :: xs + end function getTotalXS + + !! + !! Update material properties at each time step + !! First update energy using simple balance, then solve for temperature, + !! then update temperature-dependent properties + !! + !! Args: + !! tallyEnergy [in] -> Energy absorbed into material + !! printUpdate [in, optional] -> Bool, if true then will print updates to screen + !! + subroutine updateMat(self, tallyEnergy, printUpdate) + import :: mgIMCMaterial, defReal, defBool + class(mgIMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + logical(defBool), intent(in), optional :: printUpdate + end subroutine updateMat + + !! + !! Return the equilibrium radiation energy density, U_r + !! + function getEmittedRad(self) result(emittedRad) + import :: mgIMCMaterial, defReal, RNG + class(mgIMCMaterial), intent(inout) :: self + !class(RNG), intent(inout) :: rand + real(defReal) :: emittedRad + end function getEmittedRad + + !! + !! Return Fleck factor + !! + function getFleck(self) result(fleck) + import :: mgIMCMaterial, defReal + class(mgIMCMaterial), intent(in) :: self + real(defReal) :: fleck + end function getFleck + + !! + !! Get temperature of material + !! + function getTemp(self) result(T) + import :: mgIMCMaterial, defReal + class(mgIMCMaterial), intent(inout) :: self + real(defReal) :: T + end function getTemp + + !! + !! Provide material with time step size + !! + !! Args: + !! dt [in] -> time step size [s] + !! + subroutine setTimeStep(self, dt) + import :: mgIMCMaterial, defReal + class(mgIMCMaterial), intent(inout) :: self + real(defReal), intent(in) :: dt + end subroutine setTimeStep + + + end interface + + + +contains + + !! + !! Return Macroscopic XSs for the material given particle + !! + !! See IMCMaterial_inter for details + !! + subroutine getMacroXSs_byP(self, xss, p) + class(mgIMCMaterial), intent(in) :: self + type(IMCMacroXSs), intent(out) :: xss + class(particle), intent(in) :: p + character(100), parameter :: Here = 'getMacroXSs_byP (mgIMCMateerial_inter.f90)' + + if( p % isMG) then + call self % getMacroXSs(xss, p % G, p % pRNG) + + else + call fatalError(Here, 'CE particle was given to MG data') + + end if + end subroutine getMacroXSs_byP + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(mgIMCMaterial), intent(inout) :: self + + end subroutine kill + + !! + !! Cast materialHandle pointer to mgIMCMaterial pointer + !! + !! Args: + !! source [in] -> source pointer of class materialHandle + !! + !! Result: + !! Null is source is not of ceIMCMaterial + !! Pointer to source if source is ceIMCMaterial class + !! + pure function mgIMCMaterial_CptrCast(source) result(ptr) + class(materialHandle), pointer, intent(in) :: source + class(mgIMCMaterial), pointer :: ptr + + select type(source) + class is(mgIMCMaterial) + ptr => source + + class default + ptr => null() + end select + + end function mgIMCMaterial_CptrCast + +end module mgIMCMaterial_inter diff --git a/NuclearData/nuclearDataReg_mod.f90 b/NuclearData/nuclearDataReg_mod.f90 index d1f47dc72..f439218c3 100644 --- a/NuclearData/nuclearDataReg_mod.f90 +++ b/NuclearData/nuclearDataReg_mod.f90 @@ -13,6 +13,7 @@ !! Available ND TYPES: !! CE_NEUTRON !! MG_NEUTRON +!! MG_PHOTON !! !! Private members: !! databases -> Array with defined databases (name, definition, @@ -53,7 +54,7 @@ module nuclearDataReg_mod use numPrecision - use universalVariables, only : P_NEUTRON_CE, P_NEUTRON_MG + use universalVariables, only : P_NEUTRON_CE, P_NEUTRON_MG, P_PHOTON_MG use genericProcedures, only : fatalError, numToChar, printParticleType use charMap_class, only : charMap use dictionary_class, only : dictionary @@ -62,6 +63,7 @@ module nuclearDataReg_mod use nuclearDatabase_inter, only : nuclearDatabase use ceNeutronDatabase_inter, only : ceNeutronDatabase, ceNeutronDatabase_CptrCast use mgNeutronDatabase_inter, only : mgNeutronDatabase, mgNeutronDatabase_CptrCast + use mgIMCDatabase_inter, only : mgIMCDatabase, mgIMCDatabase_CptrCast use materialMenu_mod, only : mm_init => init, mm_kill => kill, mm_nMat => nMat,& mm_nameMap => nameMap @@ -74,6 +76,9 @@ module nuclearDataReg_mod ! Neutron MG use baseMgNeutronDatabase_class, only : baseMgNeutronDatabase + ! Photon MG + use baseMgIMCDatabase_class, only : baseMgIMCDatabase + implicit none private @@ -99,6 +104,7 @@ module nuclearDataReg_mod public :: kill public :: getNeutronCE public :: getNeutronMG + public :: getIMCMG public :: get public :: getMatNames @@ -113,6 +119,7 @@ module nuclearDataReg_mod character(nameLen), dimension(*), parameter :: AVAILABLE_NUCLEAR_DATABASES = & ['aceNeutronDatabase ', & 'baseMgNeutronDatabase ', & + 'baseMgIMCDatabase ', & 'aceNeutronDatabaseUni ', & 'aceNeutronDatabaseUniIdx'] @@ -126,6 +133,9 @@ module nuclearDataReg_mod class(mgNeutronDatabase), pointer :: active_mgNeutron => null() integer(shortInt) :: activeIdx_mgNeutron = 0 + class(mgIMCDatabase), pointer :: active_mgIMC => null() + integer(shortInt) :: activeIdx_mgIMC = 0 + contains !! @@ -322,6 +332,13 @@ subroutine activate(type, name, activeMat, silent) call fatalError(Here,trim(name)//' is not database for MG neutrons') end if + case(P_PHOTON_MG) + activeIdx_mgIMC = idx + active_mgIMC => mgIMCDatabase_CptrCast(ptr) + if(.not.associated(active_mgIMC)) then + call fatalError(Here,trim(name)//' is not database for MG IMC') + end if + case default call fatalError(Here,'Unrecognised type of data to activate. Check parameters. Got: '//& numToChar(type)) @@ -357,11 +374,18 @@ subroutine display() if(idx /= 0) activeName = databases(idx) % name print '(A)', " MG NEUTRON DATA: " // trim(activeName) + ! MG IMC + activename = 'NONE' + idx = activeIdx_mgIMC + if(idx /= 0) activeName = databases(idx) % name + print '(A)', " MG IMC DATA: " // trim(activeName) + ! INACTIVE DATABASES print '(A)', "INACTIVE DATABASES:" do idx=1,size(databases) if(idx == activeIdx_mgNeutron) cycle if(idx == activeIdx_ceNeutron) cycle + if(idx == activeIdx_mgIMC) cycle end do print '(A)',repeat('\/',30) @@ -397,6 +421,10 @@ subroutine kill() activeIdx_mgNeutron = 0 active_mgNeutron => null() + ! MG IMC + activeIdx_mgIMC = 0 + active_mgIMC => null() + end subroutine kill !! @@ -419,13 +447,13 @@ function getNeutronCE() result(ptr) end function getNeutronCE !! - !! Return pointer to an active Neutron CE Database + !! Return pointer to an active Neutron MG Database !! !! Args: !! None !! !! Result: - !! ceNeutronDatabase class pointer + !! mgNeutronDatabase class pointer !! !! Errors: !! If there is no active database returns NULL ptr @@ -437,6 +465,25 @@ function getNeutronMG() result(ptr) end function getNeutronMG + !! + !! Return pointer to an active IMC MG Database + !! + !! Args: + !! None + !! + !! Result: + !! mgIMCDatabase class pointer + !! + !! Errors: + !! If there is no active database returns NULL ptr + !! + function getIMCMG() result(ptr) + class(mgIMCDatabase), pointer :: ptr + + ptr => active_mgIMC + + end function getIMCMG + !! !! Return pointer to an active Nuclear Database given particle type !! @@ -463,6 +510,9 @@ function get_byType(type, where) result(ptr) case(P_NEUTRON_MG) ptr => getNeutronMG() + case(P_PHOTON_MG) + ptr => getIMCMG() + case default ptr => null() end select @@ -576,6 +626,9 @@ subroutine new_nuclearDatabase(database, type) case('baseMgNeutronDatabase') allocate(baseMgNeutronDatabase :: database) + case('baseMgIMCDatabase') + allocate(baseMgIMCDatabase :: database) + case('aceNeutronDatabaseUni') allocate(aceNeutronDatabaseUni :: database) diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index 0c123df8d..3eaff28f2 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -225,7 +225,7 @@ end function getNuclide !! if MT < 0 then reaction is associated with material: idx -> matIdx !! if MT > 0 then reaction is associated with nuclide: idx -> nucIdx !! - !! NOTE: This function can be used to enquire abou the presence of data. If the data is + !! NOTE: This function can be used to enquire about the presence of data. If the data is !! not present null() pointer is always returned! !! !! Args: diff --git a/NuclearData/xsPackages/CMakeLists.txt b/NuclearData/xsPackages/CMakeLists.txt index 0c91759d9..4993681e3 100644 --- a/NuclearData/xsPackages/CMakeLists.txt +++ b/NuclearData/xsPackages/CMakeLists.txt @@ -1 +1,2 @@ -add_sources(./neutronXsPackages_class.f90) \ No newline at end of file +add_sources(./neutronXsPackages_class.f90 + ./IMCXsPackages_class.f90) diff --git a/NuclearData/xsPackages/IMCXsPackages_class.f90 b/NuclearData/xsPackages/IMCXsPackages_class.f90 new file mode 100644 index 000000000..f124a160d --- /dev/null +++ b/NuclearData/xsPackages/IMCXsPackages_class.f90 @@ -0,0 +1,101 @@ +!! +!! This module brakes standard rules +!! It contains a library of XS Packages for IMC particle type +!! +!! Pretty much a copy of neutronXsPackages_class, may contain unnecessary lines +!! +module IMCXsPackages_class + + use numPrecision + use endfConstants + + implicit none + private + + !! + !! IMC MACROscopic Reaction XSS + !! + !! Public Members: + !! total -> total opacity [1/cm] + !! elasticScatter -> sum of MT=2 elastic IMC scattering [1/cm] + !! inelasticScatter -> sum of all IMC producing reaction that are not elastic scattering [1/cm] + !! capture -> sum of all reactions without secondary photons [1/cm] + !! planck -> frequency-normalised planck opacity of mateirial [1/cm] + !! + !! Interface: + !! clean -> Set all XSs to 0.0 + !! get -> Return XS by MT number + !! + type, public :: IMCMacroXSs + real(defReal) :: total = ZERO + real(defReal) :: elasticScatter = ZERO + real(defReal) :: inelasticScatter = ZERO + real(defReal) :: capture = ZERO + real(defReal) :: planck = ZERO + contains + procedure :: clean + procedure :: get + end type IMCMacroXSs + +contains + + !! + !! Clean IMC MacroXSs + !! + !! Sets all XSs to 0.0 + !! + !! Args: + !! None + !! + !! Errors: + !! None + !! + elemental subroutine clean(self) + class(IMCMacroXSs), intent(inout) :: self + + self % total = ZERO + self % elasticScatter = ZERO + self % inelasticScatter = ZERO + self % capture = ZERO + self % planck = ZERO + + end subroutine clean + + !! + !! Return XSs by MT number + !! + !! Args: + !! MT [in] -> Requested MT number + !! + !! Result: + !! Value of the XS + !! + !! Errors: + !! Returns 0.0 for invalid MT + !! + elemental function get(self, MT) result(xs) + class(IMCMacroXSs), intent(in) :: self + integer(shortInt), intent(in) :: MT + real(defReal) :: xs + + select case(MT) + case(macroTotal) + xs = self % total + + case(macroCapture) + xs = self % capture + + case(macroEscatter) + xs = self % elasticScatter + + case(macroPlanck) + xs = self % planck + + case default + xs = ZERO + + end select + + end function get + +end module IMCXsPackages_class diff --git a/ParticleObjects/Source/CMakeLists.txt b/ParticleObjects/Source/CMakeLists.txt index 601c905be..dcb7ce818 100644 --- a/ParticleObjects/Source/CMakeLists.txt +++ b/ParticleObjects/Source/CMakeLists.txt @@ -1,7 +1,9 @@ # Add Source Files to the global list add_sources( source_inter.f90 configSource_inter.f90 - sourceFactory_func.f90 - pointSource_class.f90 + sourceFactory_func.f90 + pointSource_class.f90 fissionSource_class.f90 + IMCSource_class.f90 + bbSurfaceSource_class.f90 ) diff --git a/ParticleObjects/Source/IMCSource_class.f90 b/ParticleObjects/Source/IMCSource_class.f90 new file mode 100644 index 000000000..0cd06b95c --- /dev/null +++ b/ParticleObjects/Source/IMCSource_class.f90 @@ -0,0 +1,352 @@ +module IMCSource_class + + use numPrecision + use endfConstants + use universalVariables + use genericProcedures, only : fatalError, rotateVector + use dictionary_class, only : dictionary + use RNG_class, only : RNG + + use particle_class, only : particle, particleState, P_PHOTON + use particleDungeon_class, only : particleDungeon + use source_inter, only : source, kill_super => kill + + use geometry_inter, only : geometry + use IMCMaterial_inter, only : IMCMaterial, IMCMaterial_CptrCast + use nuclearDataReg_mod, only : ndReg_getIMCMG => getIMCMG + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase + use materialMenu_mod, only : mm_nMat => nMat, & + mm_matName => matName + + implicit none + private + + integer(shortInt), parameter :: REJ = 1, FAST = 2 + + !! + !! IMC Source for uniform generation of photons within a material + !! + !! Angular distribution is isotropic. + !! + !! Private members: + !! isMG -> is the source multi-group? (default = .true.) + !! bottom -> Bottom corner (x_min, y_min, z_min) + !! top -> Top corner (x_max, y_max, z_max) + !! G -> Group (default = 1) + !! + !! Interface: + !! source_inter Interface + !! + !! SAMPLE INPUT: + !! imcSource { type IMCSource; } + !! + type, public,extends(source) :: imcSource + private + logical(defBool) :: isMG = .true. + real(defReal), dimension(3) :: bottom = ZERO + real(defReal), dimension(3) :: top = ZERO + real(defReal), dimension(3) :: latPitch = ZERO + integer(shortInt), dimension(3) :: latSizeN = 0 + integer(shortInt) :: G = 0 + real(defReal), dimension(6) :: bounds = ZERO + integer(shortInt) :: method = REJ + contains + procedure :: init + procedure :: append + procedure :: sampleParticle + procedure, private :: sampleIMC + procedure, private :: getMatBounds + procedure :: kill + end type imcSource + +contains + + !! + !! Initialise IMC Source + !! + !! See source_inter for details + !! + subroutine init(self, dict, geom) + class(imcSource), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), pointer, intent(in) :: geom + character(nameLen) :: method + character(100), parameter :: Here = 'init (imcSource_class.f90)' + + call dict % getOrDefault(self % G, 'G', 1) + + ! Provide geometry info to source + self % geom => geom + + ! Set bounding region + self % bounds = self % geom % bounds() + + ! Select method for position sampling + call dict % getOrDefault(method, 'method', 'rejection') + select case(method) + case('rejection') + self % method = REJ + + case('fast') + self % method = FAST + ! Get lattice dimensions + self % latSizeN = self % geom % latSizeN() + self % latPitch = (self % bounds(4:6) - self % bounds(1:3)) / self % latSizeN + + case default + call fatalError(Here, 'Unrecognised method. Should be "rejection" or "fast"') + end select + + end subroutine init + + + !! + !! Generate N particles to add to a particleDungeon without overriding + !! particles already present. + !! + !! Args: + !! dungeon [inout] -> particle dungeon to be added to + !! n [in] -> number of particles to place in dungeon + !! rand [inout] -> particle RNG object + !! + !! Result: + !! A dungeon populated with N particles sampled from the source, plus particles + !! already present in dungeon + !! + subroutine append(self, dungeon, N, rand) + class(imcSource), intent(inout) :: self + type(particleDungeon), intent(inout) :: dungeon + integer(shortInt), intent(in) :: N + class(RNG), intent(inout) :: rand + real(defReal), dimension(6) :: bounds + integer(shortInt) :: matIdx, i, Ntemp + real(defReal) :: energy, totalEnergy + type(RNG) :: pRand + class(mgIMCDatabase), pointer :: nucData + character(100), parameter :: Here = "append (IMCSource_class.f90)" + + ! Get pointer to appropriate nuclear database + nucData => ndReg_getIMCMG() + if(.not.associated(nucData)) call fatalError(Here, 'Failed to retrieve Nuclear Database') + + ! Obtain total energy + totalEnergy = nucData % getEmittedRad() + + ! Loop through materials + do matIdx = 1, mm_nMat() + + ! Get energy to be emitted from material matIdx + energy = nucData % getEmittedRad(matIdx) + + ! Choose particle numbers in proportion to material energy + if (energy > ZERO) then + Ntemp = int(N * energy / totalEnergy) + ! Enforce at least 1 particle + if (Ntemp == 0) Ntemp = 1 + + ! Set bounds for sampling + if (self % method == FAST) then + bounds = self % getMatBounds(matIdx) + else + bounds = self % bounds + end if + + ! Find energy per particle + energy = energy / Ntemp + + ! Sample particles + !$omp parallel + pRand = rand + !$omp do private(pRand) + do i=1, Ntemp + call pRand % stride(i) + call dungeon % detain(self % sampleIMC(pRand, matIdx, energy, bounds)) + end do + !$omp end do + !$omp end parallel + + end if + end do + + end subroutine append + + + !! + !! Sample particle's phase space co-ordinates + !! + !! See source_inter for details + !! + function sampleParticle(self, rand) result(p) + class(imcSource), intent(inout) :: self + class(RNG), intent(inout) :: rand + type(particleState) :: p + character(100), parameter :: Here = 'sampleParticle (IMCSource_class.f90)' + + ! Should not be called, useful to have extra inputs so use sampleIMC instead + call fatalError(Here, 'Should not be called, sampleIMC should be used instead.') + + ! Avoid compiler warning + p % G = self % G + + end function sampleParticle + + + !! + !! Sample particle's phase space co-ordinates + !! + !! Args: + !! rand [in] -> RNG + !! matIdx [in] -> index of material being sampled from + !! energy [in] -> energy-weight of sampled particle + !! bounds [in] -> bounds for position search, will be bounds of entire geometry if using + !! rejection sampling method, and bounds of single material if using fast + !! + function sampleIMC(self, rand, targetMatIdx, energy, bounds) result(p) + class(imcSource), intent(inout) :: self + class(RNG), intent(inout) :: rand + integer(shortInt), intent(in) :: targetMatIdx + real(defReal), intent(in) :: energy + real(defReal), dimension(6), intent(in) :: bounds + type(particleState) :: p + real(defReal), dimension(3) :: bottom, top, r, dir, rand3 + real(defReal) :: mu, phi + integer(shortInt) :: i, matIdx, uniqueID + character(100), parameter :: Here = 'sampleIMC (IMCSource_class.f90)' + + ! Sample particle position + bottom = bounds(1:3) + top = bounds(4:6) + i = 0 + rejection:do + rand3(1) = rand % get() + rand3(2) = rand % get() + rand3(3) = rand % get() + r = (top - bottom) * rand3 + bottom + + ! Find material under position + call self % geom % whatIsAt(matIdx, uniqueID, r) + + ! Exit if in desired material + if (matIdx == targetMatIdx) exit rejection + + ! Should exit immediately if using fast method as bounds should contain only matIdx + if (self % method == FAST) call fatalError(Here, 'Fast sourcing returned incorrect material') + + ! Protect against infinite loop + i = i+1 + if (i > 10000) call fatalError(Here, '10,000 failed attempts in rejection sampling') + + end do rejection + + ! Sample direction - chosen uniformly inside unit sphere + mu = 2 * rand % get() - 1 + phi = rand % get() * 2*pi + dir(1) = mu + dir(2) = sqrt(1-mu**2) * cos(phi) + dir(3) = sqrt(1-mu**2) * sin(phi) + + ! Assign basic phase-space coordinates + p % matIdx = matIdx + p % uniqueID = uniqueID + p % time = ZERO + p % type = P_PHOTON + p % r = r + p % dir = dir + p % G = self % G + p % isMG = .true. + p % wgt = energy + + end function sampleIMC + + + !! + !! Get location of material in lattice for position sampling + !! + !! Note that this may be incorrect depending on how lattice input is given, this function + !! assumes that geometry has been generated by discretiseGeom_class.f90 + !! + !! Args: + !! matIdx [in] -> matIdx for which to calculate bounds + !! matBounds [out] -> boundary of lattice cell, [xmin,ymin,zmin,xmax,ymax,zmax] + !! + !! TODO: + !! Would be nice to have most of this in a geometry module + !! + function getMatBounds(self, matIdx) result(matBounds) + class(imcSource), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + real(defReal), dimension(6) :: matBounds + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: i, latIdFlipped + character(nameLen) :: matName + character(100), parameter :: Here = 'getMatBounds (imcSourceClass.f90)' + + ! Extract lattice position from mat name (e.g. "m106 -> 106") + ! This is different from localID in latUniverse_class as is counting from a different + ! corner (see get_ijk function description below) + matName = mm_matName(matIdx) + read (matName(2:), '(I10)') latIdFlipped + + ! Set bounds of lattice cell containing matIdx + ijk = get_ijk(latIdFlipped, self % latSizeN) + + do i=1, 3 + matBounds(i) = (ijk(i)-1) * self % latPitch(i) + self % bounds(i) + SURF_TOL + matBounds(i+3) = ijk(i) * self % latPitch(i) + self % bounds(i) - SURF_TOL + end do + + end function getMatBounds + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(imcSource), intent(inout) :: self + + call kill_super(self) + + self % isMG = .true. + self % bounds = ZERO + self % G = 0 + + end subroutine kill + + !! + !! Generate ijk from flipped localID and shape + !! + !! Note that this is NOT the same as get_ijk in latUniverse_class. Lattice is built with first + !! map input as x_min, y_MAX, z_MAX cell, but localID begins at x_min, y_min, z_min cell. In + !! this module we want to find ijk from matIdx, which we convert to a flippedLocalID by + !! offsetting for void regions, which starts counting from the wrong corner. We therefore flip + !! ijk in the y and z directions in this function compared to instances of this function in other + !! modules. + !! + !! Args: + !! flippedlocalID [in] -> Local id of the cell between 1 and product(sizeN), + !! counting from wrong corner + !! sizeN [in] -> Number of cells in each cardinal direction x, y & z + !! + !! Result: + !! Array ijk which has integer position in each cardinal direction + !! + pure function get_ijk(flippedLocalID, sizeN) result(ijk) + integer(shortInt), intent(in) :: flippedLocalID + integer(shortInt), dimension(3), intent(in) :: sizeN + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: temp, base + + temp = flippedLocalID - 1 + base = temp / sizeN(1) + ijk(1) = temp - sizeN(1) * base + 1 + + temp = base + base = temp / sizeN(2) + ijk(2) = sizeN(2)*(1 + base) - temp + + ijk(3) = sizeN(3) - base + + end function get_ijk + + +end module IMCSource_class diff --git a/ParticleObjects/Source/bbSurfaceSource_class.f90 b/ParticleObjects/Source/bbSurfaceSource_class.f90 new file mode 100644 index 000000000..81f116710 --- /dev/null +++ b/ParticleObjects/Source/bbSurfaceSource_class.f90 @@ -0,0 +1,297 @@ +module bbSurfaceSource_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, numToChar + use particle_class, only : particleState, P_NEUTRON, P_PHOTON + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary + use configSource_inter, only : configSource, kill_super => kill + use geometry_inter, only : geometry + use RNG_class, only : RNG + + implicit none + private + + !! + !! Generates a source representing a black body surface + !! + !! Private members: + !! r -> bottom corner of source + !! dr -> size of surface, will be 0 in one dimension + !! dir -> direction of dominant movement: [1,0,0], [-1,0,0], [0,1,0], etc. + !! particleType -> source particle type (photon) + !! isMG -> is the source multi-group? (yes) + !! + !! Interface: + !! init -> initialise point source + !! append -> source particles and add to existing dungeon + !! sampleType -> set particle type + !! samplePosition -> set particle position + !! sampleEnergyAngle -> sample particle angle + !! sampleEnergy -> set particle energy (isMG = .true., G = 1) + !! sampleWeight -> set particle energy-weight + !! kill -> terminate source + !! + !! Sample Dictionary Input: + !! source { + !! type bbSurfaceSource; + !! r (x_min x_max y_min y_max z_min z_max); -> Position bounds of surface + !! -> min and max must be equal in one dimension + !! #dir -1; -> optional, negative will reverse direction in dominant axis + !! -> defaults to positive + !! temp 1; -> temperature of the black body source + !! #deltaT 0.05; -> time step size, automatically added to dictionary in IMCPhysicsPackage_class.f90 + !! N 100; -> number of particles per time step, only used if append is called with N = 0 + !! } + !! + type, public,extends(configSource) :: bbSurfaceSource + private + real(defReal), dimension(3) :: r = ZERO + real(defReal), dimension(3) :: dr = ZERO + integer(shortInt), dimension(3) :: dir = ZERO + integer(shortInt) :: particleType = P_PHOTON + logical(defBool) :: isMG = .true. + real(defReal) :: T = ZERO + real(defReal) :: deltaT = ZERO + integer(shortInt) :: N = 0 + contains + procedure :: init + procedure :: append + procedure :: sampleType + procedure :: samplePosition + procedure :: sampleEnergy + procedure :: sampleEnergyAngle + procedure :: sampleWeight + procedure :: kill + end type bbSurfaceSource + +contains + + !! + !! Initialise from dictionary + !! + !! See source_inter for details + !! + !! Errors: + !! - error if an axis other than x, y, or z is given + !! + subroutine init(self, dict, geom) + class(bbSurfaceSource), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), pointer, intent(in) :: geom + real(defReal), dimension(:), allocatable :: temp + integer(shortInt) :: i, dir + character(100), parameter :: Here = 'init (bbSurfaceSource_class.f90)' + + ! Provide geometry info to source + self % geom => geom + + ! Provide particle type + self % particleType = P_PHOTON + + ! Get and check position vector + call dict % get(temp, 'r') + if (size(temp) /= 6) call fatalError(Here, 'r should be of size 6') + do i = 1, 3 + ! Store x_min, y_min, z_min + self % r(i) = temp(2*i-1) + ! Store dx, dy, dz + self % dr(i) = temp(2*i) - temp(2*i-1) + ! Check for compatible min and max + if (self % dr(i) < 0) call fatalError(Here, 'Min > Max along direction '//numToChar(i)) + end do + ! Check that exactly one normal axis is present + if (count(self % dr == 0) /= 1) call fatalError(Here, 'No clearly defined axis extracted') + + ! Get primary direction + call dict % getOrDefault(dir, 'dir', 1) + do i = 1, 3 + if (self % dr(i) == 0) self % dir(i) = sign(1, dir) + end do + + ! Move by 2*SURF_TOL to ensure sourcing in correct material + self % r = self % r + 2*SURF_TOL*self % dir + + ! Get remaining information + call dict % get(self % T, 'temp') + call dict % get(self % deltaT, 'deltaT') ! Automatically added to dict in IMC physics package + call dict % getOrDefault(self % N, 'N', 1) + + end subroutine init + + !! + !! Add particles to given dungeon + !! + !! See source_inter for details + !! + !! If N is given as 0, then N is instead taken from the input dictionary defining this source + !! to allow PP to have control over particle numbers + !! + subroutine append(self, dungeon, N, rand) + class(bbSurfaceSource), intent(inout) :: self + type(particleDungeon), intent(inout) :: dungeon + integer(shortInt), intent(in) :: N + class(RNG), intent(inout) :: rand + integer(shortInt) :: i + type(RNG) :: pRand + character(100), parameter :: Here = 'append (bbSurfaceSource_class.f90)' + + ! Set number to generate. Using 0 in function call will use N from input dictionary + if (N /= 0) self % N = N + ! TODO change so that this override is only temporary, so that can be called with 0 again later + + ! Generate N particles to populate dungeon + !$omp parallel + pRand = rand + !$omp do private(pRand) + do i = 1, self % N + call pRand % stride(i) + call dungeon % detain(self % sampleParticle(pRand)) + end do + !$omp end do + !$omp end parallel + + end subroutine append + + !! + !! Provide particle type + !! + !! See configSource_inter for details. + !! + subroutine sampleType(self, p, rand) + class(bbSurfaceSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + p % type = self % particleType + + end subroutine sampleType + + !! + !! Provide particle position + !! + !! See configSource_inter for details. + !! + subroutine samplePosition(self, p, rand) + class(bbSurfaceSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + integer(shortInt) :: i + real(defReal), dimension(3) :: r + + ! Set new x, y and z coords + do i = 1, 3 + r(i) = (self % dr(i)) * rand % get() + self % r(i) + end do + + ! Assign to particle + p % r = r + + end subroutine samplePosition + + !! + !! Sample angle + !! + !! See configSource_inter for details. + !! + subroutine sampleEnergyAngle(self, p, rand) + class(bbSurfaceSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + real(defReal), dimension(3) :: dir + real(defReal) :: phi, mu + character(100), parameter :: Here = 'sampleEnergyAngle (bbSurfaceSource_class.f90)' + + ! Sample required phi and mu + phi = TWO_PI * rand % get() + mu = sqrt(rand % get()) + + ! Choose direction based on dominant direction given in self % dir + if (self % dir(1) == 1) then ! Positive x + dir = [ mu, sqrt(1-mu*mu)*cos(phi), sqrt(1-mu*mu)*sin(phi)] + + else if (self % dir(1) == -1) then ! Negative x + dir = [-mu, sqrt(1-mu*mu)*cos(phi), sqrt(1-mu*mu)*sin(phi)] + + else if (self % dir(2) == 1) then ! Positive y + dir = [sqrt(1-mu*mu)*sin(phi), mu, sqrt(1-mu*mu)*cos(phi)] + + else if (self % dir(2) == -1) then ! Negative y + dir = [sqrt(1-mu*mu)*sin(phi), -mu, sqrt(1-mu*mu)*cos(phi)] + + else if (self % dir(3) == 1) then ! Positive z + dir = [sqrt(1-mu*mu)*cos(phi), sqrt(1-mu*mu)*sin(phi), mu] + + else if (self % dir(3) == -1) then ! Negative z + dir = [sqrt(1-mu*mu)*cos(phi), sqrt(1-mu*mu)*sin(phi), -mu] + + else + call fatalError(Here, 'Invalid direction vector') + end if + + ! Assign to particle + p % dir = dir + + end subroutine sampleEnergyAngle + + !! + !! Provide particle energy, currently only a single group + !! + !! See configSource_inter for details. + !! + subroutine sampleEnergy(self, p, rand) + class(bbSurfaceSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + p % isMG = .true. + p % G = 1 + + end subroutine sampleEnergy + + !! + !! Provide particle energy-weight + !! + !! Sampled as a black body surface, see "Four Decades of Implicit Monte Carlo", + !! Allan B Wollaber, p.24-25 + !! + !! See configSource_inter for details. + !! + subroutine sampleWeight(self, p, rand) + class(bbSurfaceSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + real(defReal) :: area, num + + ! Calculate surface area of source + area = product(self % dr, self % dr /= ZERO) + + ! Calculate energy weight per particle + num = radiationConstant * lightSpeed * self % deltaT * self % T**4 * area + p % wgt = num / (4 * self % N) + + end subroutine sampleWeight + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(bbSurfaceSource), intent(inout) :: self + + ! Kill superclass + call kill_super(self) + + ! Kill local components + self % r = ZERO + self % dr = ZERO + self % dir = ZERO + self % particleType = P_PHOTON + self % isMG = .true. + self % T = ZERO + self % deltaT = ZERO + self % N = ZERO + + end subroutine kill + +end module bbSurfaceSource_class diff --git a/ParticleObjects/Source/configSource_inter.f90 b/ParticleObjects/Source/configSource_inter.f90 index 968dd0dbe..6aa274c80 100644 --- a/ParticleObjects/Source/configSource_inter.f90 +++ b/ParticleObjects/Source/configSource_inter.f90 @@ -32,6 +32,7 @@ module configSource_inter contains procedure :: sampleParticle + procedure :: sampleWeight procedure(sampleType), deferred :: sampleType procedure(samplePosition), deferred :: samplePosition procedure(sampleEnergy), deferred :: sampleEnergy @@ -133,11 +134,24 @@ function sampleParticle(self, rand) result(p) call self % samplePosition(p, rand) call self % sampleEnergyAngle(p, rand) call self % sampleEnergy(p, rand) + call self % sampleWeight(p, rand) p % time = ZERO - p % wgt = ONE end function sampleParticle + !! + !! Set particle's weight to 1. + !! Can be overriden in subclasses if needed. + !! + subroutine sampleWeight(self, p, rand) + class(configSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + p % wgt = ONE + + end subroutine sampleWeight + !! !! Return to uninitialised state !! diff --git a/ParticleObjects/Source/sourceFactory_func.f90 b/ParticleObjects/Source/sourceFactory_func.f90 index b4e15c429..91c7ec6d3 100644 --- a/ParticleObjects/Source/sourceFactory_func.f90 +++ b/ParticleObjects/Source/sourceFactory_func.f90 @@ -8,8 +8,10 @@ module sourceFactory_func use source_inter, only : source ! source implementations - use pointSource_class, only : pointSource - use fissionSource_class, only : fissionSource + use pointSource_class, only : pointSource + use fissionSource_class, only : fissionSource + use IMCSource_class, only : imcSource + use bbSurfaceSource_class, only : bbSurfaceSource ! geometry use geometry_inter, only : geometry @@ -24,8 +26,10 @@ module sourceFactory_func ! It is printed if type was unrecognised ! NOTE: ! For now it is necessary to adjust trailing blanks so all entries have the same length - character(nameLen),dimension(*),parameter :: AVAILABLE_sources = [ 'pointSource ',& - 'fissionSource'] + character(nameLen),dimension(*),parameter :: AVAILABLE_sources = [ 'pointSource ',& + 'fissionSource ',& + 'imcSource ',& + 'bbSurfaceSource'] contains @@ -57,6 +61,14 @@ subroutine new_source(new, dict, geom) allocate(fissionSource :: new) call new % init(dict, geom) + case('imcSource') + allocate(imcSource :: new) + call new % init(dict, geom) + + case('bbSurfaceSource') + allocate(bbSurfaceSource :: new) + call new % init(dict, geom) + !*** NEW SOURCE TEMPLATE ***! !case('') ! allocate( :: new) diff --git a/ParticleObjects/Source/source_inter.f90 b/ParticleObjects/Source/source_inter.f90 index a67ba9c63..ae3254059 100644 --- a/ParticleObjects/Source/source_inter.f90 +++ b/ParticleObjects/Source/source_inter.f90 @@ -1,11 +1,12 @@ module source_inter use numPrecision - use particle_class, only : particleState + use particle_class, only : particle, particleState use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary use RNG_class, only : RNG use geometry_inter, only : geometry + use genericProcedures, only : fatalError implicit none private @@ -29,6 +30,7 @@ module source_inter !! Interface: !! init -> initialise the source !! generate -> generate particles to fill a dungeon + !! append -> generate new particles to add to an existing dungeon !! sampleParticle -> sample particles from the corresponding distributions !! kill -> clean up the source !! @@ -36,7 +38,8 @@ module source_inter private class(geometry), pointer, public :: geom => null() contains - procedure, non_overridable :: generate + procedure :: generate + procedure :: append procedure(sampleParticle), deferred :: sampleParticle procedure(init), deferred :: init procedure(kill), deferred :: kill @@ -121,6 +124,36 @@ subroutine generate(self, dungeon, n, rand) end subroutine generate + !! + !! Generate particles to add to a particleDungeon without overriding + !! particles already present + !! + !! Adds to a particle dungeon N particles, sampled + !! from the corresponding source distributions + !! + !! Args: + !! dungeon [inout] -> particle dungeon to be added to + !! n [in] -> number of particles to place in dungeon + !! rand [inout] -> particle RNG object + !! + !! Result: + !! A dungeon populated with n particles sampled from the source, plus + !! particles already present in dungeon + !! + subroutine append(self, dungeon, N, rand) + class(source), intent(inout) :: self + type(particleDungeon), intent(inout) :: dungeon + integer(shortInt), intent(in) :: N + class(RNG), intent(inout) :: rand + integer(shortInt) :: i + + ! Generate n particles to populate dungeon + do i = 1, N + call dungeon % detain(self % sampleParticle(rand)) + end do + + end subroutine append + !! !! Return to uninitialised state !! diff --git a/ParticleObjects/particleDungeon_class.f90 b/ParticleObjects/particleDungeon_class.f90 index c4d532f7d..8e86a8ba5 100644 --- a/ParticleObjects/particleDungeon_class.f90 +++ b/ParticleObjects/particleDungeon_class.f90 @@ -53,6 +53,9 @@ module particleDungeon_class !! does not take nonuniform weight of particles into account !! setSize(n) -> sizes dungeon to have n dummy particles for ease of overwriting !! printToFile(name) -> prints population in ASCII format to file "name" + !! printToScreen(prop,nMax,total) -> prints property to screen for up to nMax particles + !! popSize() -> returns number of particles in dungeon + !! popWeight() -> returns total population weight !! !! Build procedures: !! init(maxSize) -> allocate space to store maximum of maxSize particles @@ -60,8 +63,8 @@ module particleDungeon_class !! type, public :: particleDungeon private - real(defReal),public :: k_eff = ONE ! k-eff for fission site generation rate normalisation - integer(shortInt) :: pop = 0 ! Current population size of the dungeon + real(defReal),public :: k_eff = ONE ! k-eff for fission site generation rate normalisation + integer(shortInt) :: pop = 0 ! Current population size of the dungeon ! Storage space type(particleState), dimension(:), allocatable :: prisoners @@ -91,6 +94,7 @@ module particleDungeon_class procedure :: popWeight procedure :: setSize procedure :: printToFile + procedure :: printToScreen ! Private procedures procedure, private :: detain_particle @@ -447,7 +451,7 @@ pure subroutine cleanPop(self) end subroutine cleanPop !! - !! Returns number of neutrons in the dungeon + !! Returns number of particles in the dungeon !! function popSize(self) result(pop) class(particleDungeon), intent(in) :: self @@ -517,14 +521,15 @@ subroutine printToFile(self, name) integer(shortInt) :: i filename = trim(name)//'.txt' - open(unit = 10, file = filename, status = 'new') + open(unit = 10, file = filename) ! Print out each particle co-ordinate do i = 1, self % pop write(10,'(8A)') numToChar(self % prisoners(i) % r), & numToChar(self % prisoners(i) % dir), & numToChar(self % prisoners(i) % E), & - numToChar(self % prisoners(i) % G) + numToChar(self % prisoners(i) % G), & + numToChar(self % prisoners(i) % matIdx) end do ! Close the file @@ -532,4 +537,93 @@ subroutine printToFile(self, name) end subroutine printToFile + !! + !! Prints given property of particles to screen + !! + !! Args: + !! prop [in] -> Particle property to be displayed + !! nMax [in] -> Maximum number of particles displayed + !! total [in] -> Optional, if True then sum contributions of particles + !! and print for total + !! + !! Errors: + !! fatalError if prop is invalid + !! + subroutine printToScreen(self, prop, nMax) + class(particleDungeon), intent(in) :: self + character(*), intent(in) :: prop + integer(shortInt), intent(in) :: nMax + integer(shortInt) :: i,iMax + character(100), parameter :: Here = 'printToScreen (particleDungeon_class.f90)' + + character(nameLen), dimension(*), parameter :: AVAILABLE_props = [ 'r ',& + 'dir ',& + 'matIdx',& + 'E ',& + 'G ',& + 'wgt ',& + 'time ',& + 'pop '] + + print *, 'Number in dungeon =', self % pop + + ! Number of particles to be printed + iMax = min(nMax, self % pop) + + ! Print desired quantities + select case(prop) + case('r') + print *, '** ** Position ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % r) + end do + + case('dir') + print *, '** ** Direction ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % dir) + end do + + case('matIdx') + print *, '** ** matIdx ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % matIdx) + end do + + case('E') + print *, '** ** Energy ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % E) + end do + + case('G') + print *, '** ** Group ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % G) + end do + + case('wgt') + print *, '** ** Weight ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % wgt) + end do + + case('time') + print *, '** ** Time ** **' + do i = 1, iMax + print *, i,numToChar(self % prisoners(i) % time) + end do + + case('pop') + ! Do nothing, pop already printed above + + case default + print *, AVAILABLE_props + call fatalError(Here, 'Unrecognised particle property : ' // trim(prop)) + + end select + + end subroutine printToScreen + + end module particleDungeon_class diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index e88bd73d6..52e252864 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -1,5 +1,6 @@ module particle_class + use numPrecision use universalVariables use genericProcedures @@ -135,6 +136,7 @@ module particle_class procedure :: getUniIdx procedure :: matIdx procedure, non_overridable :: getType + procedure :: getSpeed ! Operations on coordinates procedure :: moveGlobal @@ -398,13 +400,13 @@ pure function matIdx(self) result(Idx) end function matIdx !! - !! Return one of the particle Tpes defined in universal variables + !! Return one of the particle types defined in universal variables !! !! Args: !! None !! !! Result: - !! P_NEUTRON_CE, P_NEUTRON_MG + !! P_NEUTRON_CE, P_NEUTRON_MG, P_PHOTON_MG !! !! Errors: !! None @@ -413,14 +415,46 @@ pure function getType(self) result(type) class(particle), intent(in) :: self integer(shortInt) :: type - if (self % isMG) then - type = P_NEUTRON_MG - else - type = P_NEUTRON_CE + ! Check for neutron + if (self % type == P_NEUTRON) then + if (self % isMg) then + type = P_NEUTRON_MG + else + type = P_NEUTRON_CE + end if + ! Check for photon + else if (self % type == P_PHOTON) then + if (self % isMg) then + type = P_PHOTON_MG + else + type = P_PHOTON_CE + end if end if end function getType + !! + !! Return speed of particle + !! + !! Args: + !! None + !! + !! Errors: + !! Currently returns lightSpeed for P_PHOTONs and gives error otherwise + !! + function getSpeed(self) result(speed) + class(particle), intent(in) :: self + real(defReal) :: speed + character(100), parameter :: Here = 'getSpeed (particle_class.f90)' + + if (self % type == P_PHOTON) then + speed = lightSpeed + else + call fatalError(Here, "Not yet coded to provide speed for neutrons") + end if + + end function getSpeed + !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> !! Particle operations on coordinates procedures !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> diff --git a/PhysicsPackages/CMakeLists.txt b/PhysicsPackages/CMakeLists.txt index a1db7a1f4..c7c10a675 100644 --- a/PhysicsPackages/CMakeLists.txt +++ b/PhysicsPackages/CMakeLists.txt @@ -3,7 +3,8 @@ add_sources( ./physicsPackage_inter.f90 ./physicsPackageFactory_func.f90 ./eigenPhysicsPackage_class.f90 - ./fixedSourcePhysicsPackage_class.f90 + ./fixedSourcePhysicsPackage_class.f90 + ./IMCPhysicsPackage_class.f90 ./vizPhysicsPackage_class.f90 ./rayVolPhysicsPackage_class.f90 ) diff --git a/PhysicsPackages/IMCPhysicsPackage_class.f90 b/PhysicsPackages/IMCPhysicsPackage_class.f90 new file mode 100644 index 000000000..bd68ebd08 --- /dev/null +++ b/PhysicsPackages/IMCPhysicsPackage_class.f90 @@ -0,0 +1,541 @@ +module IMCPhysicsPackage_class + + use numPrecision + use universalVariables + use endfConstants + use genericProcedures, only : fatalError, printFishLineR, numToChar, rotateVector + use hashFunctions_func, only : FNV_1 + use dictionary_class, only : dictionary + use outputFile_class, only : outputFile + + ! Timers + use timer_mod, only : registerTimer, timerStart, timerStop, & + timerTime, timerReset, secToChar + + ! Particle classes and Random number generator + use particle_class, only : particle, P_PHOTON + use particleDungeon_class, only : particleDungeon + use source_inter, only : source + use RNG_class, only : RNG + + ! Physics package interface + use physicsPackage_inter, only : physicsPackage + + ! Geometry + use geometry_inter, only : geometry + use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_addGeom => addGeom, & + gr_geomIdx => geomIdx + use discretiseGeom_class, only : discretise + + ! Nuclear Data + use materialMenu_mod, only : mm_nMat => nMat ,& + mm_matName => matName + use nuclearDataReg_mod, only : ndReg_init => init ,& + ndReg_activate => activate ,& + ndReg_display => display, & + ndReg_kill => kill, & + ndReg_get => get ,& + ndReg_getMatNames => getMatNames + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase, mgIMCDatabase_CptrCast + use IMCMaterial_inter, only : IMCMaterial, IMCMaterial_CptrCast + + ! Operators + use collisionOperator_class, only : collisionOperator + use transportOperator_inter, only : transportOperator + + ! Tallies + use tallyCodes + use tallyAdmin_class, only : tallyAdmin + use tallyResult_class, only : tallyResult + use absorptionClerk_class, only : absClerkResult + + ! Factories + use transportOperatorFactory_func, only : new_transportOperator + use sourceFactory_func, only : new_source + + implicit none + + private + + !! + !! Physics Package for IMC calculations + !! + type, public,extends(physicsPackage) :: IMCPhysicsPackage + private + ! Building blocks +! class(nuclearDatabase), pointer :: nucData => null() + class(mgIMCDatabase), pointer :: nucData => null() + class(geometry), pointer :: geom => null() + integer(shortInt) :: geomIdx = 0 + type(collisionOperator) :: collOp + class(transportOperator), allocatable :: transOp + class(RNG), pointer :: pRNG => null() + type(tallyAdmin),pointer :: tally => null() + type(tallyAdmin),pointer :: imcWeightAtch => null() + + ! Settings + integer(shortInt) :: N_steps + integer(shortInt) :: pop + integer(shortInt) :: limit + real(defReal) :: deltaT + character(pathLen) :: outputFile + character(nameLen) :: outputFormat + integer(shortInt) :: printSource = 0 + integer(shortInt) :: particleType + logical(defBool) :: sourceGiven = .false. + integer(shortInt) :: nMat + integer(shortInt) :: printUpdates + + ! Calculation components + type(particleDungeon), pointer :: thisStep => null() + type(particleDungeon), pointer :: nextStep => null() + type(particleDungeon), pointer :: temp_dungeon => null() + class(source), allocatable :: inputSource + class(source), allocatable :: IMCSource + + ! Timer bins + integer(shortInt) :: timerMain + real (defReal) :: CPU_time_start + real (defReal) :: CPU_time_end + + contains + procedure :: init + procedure :: printSettings + procedure :: steps + procedure :: collectResults + procedure :: run + procedure :: kill + + end type IMCPhysicsPackage + +contains + + subroutine run(self) + class(IMCPhysicsPackage), intent(inout) :: self + + print *, repeat("<>",50) + print *, "/\/\ IMC CALCULATION /\/\" + + call self % steps(self % tally, self % imcWeightAtch, self % N_steps) + call self % collectResults() + + print * + print *, "\/\/ END OF IMC CALCULATION \/\/" + print * + end subroutine + + !! + !! Run steps for calculation + !! + subroutine steps(self, tally, tallyAtch, N_steps) + class(IMCPhysicsPackage), intent(inout) :: self + type(tallyAdmin), pointer,intent(inout) :: tally + type(tallyAdmin), pointer,intent(inout) :: tallyAtch + integer(shortInt), intent(in) :: N_steps + integer(shortInt) :: i, j, N, num, nParticles + type(particle), save :: p + real(defReal) :: elapsed_T, end_T, T_toEnd + real(defReal), dimension(:), allocatable :: tallyEnergy + class(IMCMaterial), pointer :: mat + character(100),parameter :: Here ='steps (IMCPhysicsPackage_class.f90)' + class(tallyResult), allocatable :: tallyRes + type(collisionOperator), save :: collOp + class(transportOperator), allocatable, save :: transOp + type(RNG), target, save :: pRNG + !$omp threadprivate(p, collOp, transOp, pRNG) + + !$omp parallel + p % geomIdx = self % geomIdx + + ! Create a collision + transport operator which can be made thread private + collOp = self % collOp + transOp = self % transOp + + !$omp end parallel + + ! Reset and start timer + call timerReset(self % timerMain) + call timerStart(self % timerMain) + + allocate(tallyEnergy(self % nMat)) + + do i=1,N_steps + + ! Update tracking grid if needed by transport operator + if (associated(self % transOp % grid)) call self % transOp % grid % update() + + ! Swap dungeons to store photons remaining from previous time step + self % temp_dungeon => self % nextStep + self % nextStep => self % thisStep + self % thisStep => self % temp_dungeon + call self % nextStep % cleanPop() + + ! Select total number of particles to generate from material emission + N = self % pop + if (N + self % thisStep % popSize() > self % limit) then + ! Fleck and Cummings IMC Paper, eqn 4.11 + N = self % limit - self % thisStep % popSize() - self % nMat - 1 + end if + + ! Add to dungeon particles emitted from material + call self % IMCSource % append(self % thisStep, N, self % pRNG) + + ! Generate from input source + if( self % sourceGiven ) then + call self % inputSource % append(self % thisStep, 0, self % pRNG) + end if + + if(self % printSource == 1) then + call self % thisStep % printToFile(trim(self % outputFile)//'_source'//numToChar(i)) + end if + + call tally % reportCycleStart(self % thisStep) + + nParticles = self % thisStep % popSize() + + !$omp parallel do schedule(dynamic) + gen: do num = 1, nParticles + + ! Create RNG which can be thread private + pRNG = self % pRNG + p % pRNG => pRNG + call p % pRNG % stride(num) + + ! Obtain paticle from dungeon + call self % thisStep % release(p) + call self % geom % placeCoord(p % coords) + + ! Check particle type + if (p % getType() /= P_PHOTON_MG) then + call fatalError(Here, 'Particle is not of type P_PHOTON_MG') + end if + + ! Assign maximum particle time + p % timeMax = self % deltaT * i + + ! For newly sourced particles, sample time uniformly within time step + if (p % time == ZERO) then + p % time = (p % pRNG % get() + i-1) * self % deltaT + end if + + ! Check for time errors + if (p % time >= p % timeMax .or. p % time < self % deltaT*(i-1)) then + call fatalError(Here, 'Particle time is not within timestep bounds') + else if (p % time /= p % time) then + call fatalError(Here, 'Particle time is NaN') + end if + + ! Save state + call p % savePreHistory() + + ! Transport particle until its death + history: do + call transOp % transport(p, tally, self % thisStep, self % nextStep) + if(p % isDead) exit history + + if(p % fate == AGED_FATE) then + ! Store particle for use in next time step + p % fate = 0 + call self % nextStep % detain(p) + exit history + end if + + call collOp % collide(p, tally, self % thisStep, self % nextStep) + + if(p % isDead) exit history + + end do history + + end do gen + !$omp end parallel do + + ! Update RNG + call self % pRNG % stride(nParticles) + + ! Send end of time step report + call tally % reportCycleEnd(self % thisStep) + + ! Calculate times + call timerStop(self % timerMain) + elapsed_T = timerTime(self % timerMain) + + ! Predict time to end + end_T = real(N_steps,defReal) * elapsed_T / i + T_toEnd = max(ZERO, end_T - elapsed_T) + + ! Display progress + call printFishLineR(i) + print * + print * + print *, 'Source batch: ', numToChar(i), ' of ', numToChar(N_steps) + print *, 'Pop: ', numToChar(self % nextStep % popSize()) + print *, 'Elapsed time: ', trim(secToChar(elapsed_T)) + print *, 'End time: ', trim(secToChar(end_T)) + print *, 'Time to end: ', trim(secToChar(T_toEnd)) + call tally % display() + + ! Obtain energy deposition tally results + call tallyAtch % getResult(tallyRes, 'imcWeightTally') + + select type(tallyRes) + class is(absClerkResult) + do j = 1, self % nMat + tallyEnergy(j) = tallyRes % clerkResults(j) + end do + class default + call fatalError(Here, 'Tally result class should be absClerkResult') + end select + + ! Update material properties + call self % nucData % updateProperties(tallyEnergy, self % printUpdates) + + ! Reset tally for next time step + call tallyAtch % reset('imcWeightTally') + + end do + + ! Output final mat temperatures + open(unit = 10, file = 'temps.txt') + do j = 1, self % nMat + mat => IMCMaterial_CptrCast(self % nucData % getMaterial(j)) + write(10, '(8A)') mm_matName(j), numToChar(mat % getTemp()) + end do + close(10) + + end subroutine steps + + !! + !! Print calculation results to file + !! + subroutine collectResults(self) + class(IMCPhysicsPackage), intent(inout) :: self + type(outputFile) :: out + character(nameLen) :: name + + call out % init(self % outputFormat) + + name = 'seed' + call out % printValue(self % pRNG % getSeed(),name) + + name = 'pop' + call out % printValue(self % pop,name) + + name = 'Source_batches' + call out % printValue(self % N_steps,name) + + call cpu_time(self % CPU_time_end) + name = 'Total_CPU_Time' + call out % printValue((self % CPU_time_end - self % CPU_time_start),name) + + name = 'Transport_time' + call out % printValue(timerTime(self % timerMain),name) + + ! Print tally + call self % tally % print(out) + + call out % writeToFile(self % outputFile) + + end subroutine collectResults + + + !! + !! Initialise from individual components and dictionaries for source and tally + !! + subroutine init(self, dict) + class(IMCPhysicsPackage), intent(inout) :: self + class(dictionary), intent(inout) :: dict + class(dictionary), pointer :: tempDict, geomDict, dataDict + type(dictionary) :: locDict1, locDict2, locDict3, locDict4 + integer(shortInt) :: seed_temp + integer(longInt) :: seed + character(10) :: time + character(8) :: date + character(:),allocatable :: string + character(nameLen) :: nucData, geomName + type(outputFile) :: test_out + integer(shortInt) :: i + character(nameLen), dimension(:), allocatable :: mats + integer(shortInt), dimension(:), allocatable :: latSizeN + type(dictionary),target :: newGeom, newData + character(100), parameter :: Here ='init (IMCPhysicsPackage_class.f90)' + + call cpu_time(self % CPU_time_start) + + ! Read calculation settings + call dict % get(self % pop,'pop') + call dict % get(self % limit, 'limit') + call dict % get(self % N_steps,'steps') + call dict % get(self % deltaT,'timeStepSize') + call dict % getOrDefault(self % printUpdates, 'printUpdates', 0) + self % particleType = P_PHOTON_MG + nucData = 'mg' + + ! Read outputfile path + call dict % getOrDefault(self % outputFile,'outputFile','./output') + + ! Get output format and verify + ! Initialise output file before calculation (so mistake in format will be cought early) + call dict % getOrDefault(self % outputFormat, 'outputFormat', 'asciiMATLAB') + call test_out % init(self % outputFormat) + + ! Register timer + self % timerMain = registerTimer('transportTime') + + ! Initialise RNG + allocate(self % pRNG) + + ! *** It is a bit silly but dictionary cannot store longInt for now + ! so seeds are limited to 32 bits (can be -ve) + if( dict % isPresent('seed')) then + call dict % get(seed_temp,'seed') + + else + ! Obtain time string and hash it to obtain random seed + call date_and_time(date, time) + string = date // time + call FNV_1(string,seed_temp) + + end if + seed = seed_temp + call self % pRNG % init(seed) + + ! Read whether to print particle source each time step + call dict % getOrDefault(self % printSource, 'printSource', 0) + + ! Automatically split geometry into a uniform grid + if (dict % isPresent('discretise')) then + + ! Store dimensions of lattice + tempDict => dict % getDictPtr('discretise') + call tempDict % get(latSizeN, 'dimensions') + + ! Create new input + call discretise(dict, newGeom, newData) + + geomDict => newGeom + dataDict => newData + + else + geomDict => dict % getDictPtr("geometry") + dataDict => dict % getDictPtr("nuclearData") + + end if + + ! Build Nuclear Data + call ndReg_init(dataDict) + + ! Build geometry + geomName = 'IMCGeom' + call gr_addGeom(geomName, geomDict) + self % geomIdx = gr_geomIdx(geomName) + self % geom => gr_geomPtr(self % geomIdx) + + ! Activate Nuclear Data *** All materials are active + call ndReg_activate(self % particleType, nucData, self % geom % activeMats()) + self % nucData => mgIMCDatabase_CptrCast(ndReg_get(self % particleType)) + + call newGeom % kill() + call newData % kill() + + ! Initialise IMC source + if (dict % isPresent('matSource')) then + tempDict => dict % getDictPtr('matSource') + call new_source(self % IMCSource, tempDict, self % geom) + else + call locDict1 % init(1) + call locDict1 % store('type', 'imcSource') + call new_source(self % IMCSource, locDict1, self % geom) + call locDict1 % kill() + end if + + ! Read external particle source definition + if( dict % isPresent('source') ) then + tempDict => dict % getDictPtr('source') + call tempDict % store('deltaT', self % deltaT) + call new_source(self % inputSource, tempDict, self % geom) + self % sourceGiven = .true. + end if + + ! Build collision operator + tempDict => dict % getDictPtr('collisionOperator') + call self % collOp % init(tempDict) + + ! Build transport operator + tempDict => dict % getDictPtr('transportOperator') + call new_transportOperator(self % transOp, tempDict) + + ! Initialise tally Admin + tempDict => dict % getDictPtr('tally') + allocate(self % tally) + call self % tally % init(tempDict) + + ! Provide materials with time step + call self % nucData % setTimeStep(self % deltaT) + + ! Store number of materials + self % nMat = mm_nMat() + self % printUpdates = min(self % printUpdates, self % nMat) + + ! Create array of material names + allocate(mats(self % nMat)) + do i=1, self % nMat + mats(i) = mm_matName(i) + end do + + ! Initialise imcWeight tally attachment + call locDict1 % init(1) + call locDict2 % init(4) + call locDict3 % init(2) + call locDict4 % init(1) + + call locDict4 % store('type', 'weightResponse') + call locDict3 % store('type','materialMap') + call locDict3 % store('materials', [mats]) + call locDict2 % store('response', ['imcWeightResponse']) + call locDict2 % store('imcWeightResponse', locDict4) + call locDict2 % store('type','absorptionClerk') + call locDict2 % store('map', locDict3) + call locDict1 % store('imcWeightTally', locDict2) + + allocate(self % imcWeightAtch) + call self % imcWeightAtch % init(locDict1) + + call self % tally % push(self % imcWeightAtch) + + ! Size particle dungeons + allocate(self % thisStep) + call self % thisStep % init(self % limit) + allocate(self % nextStep) + call self % nextStep % init(self % limit) + + call self % printSettings() + + end subroutine init + + !! + !! Deallocate memory + !! + subroutine kill(self) + class(IMCPhysicsPackage), intent(inout) :: self + + ! TODO: This subroutine + + end subroutine kill + + !! + !! Print settings of the physics package + !! + subroutine printSettings(self) + class(IMCPhysicsPackage), intent(in) :: self + + print *, repeat("<>",50) + print *, "/\/\ IMC CALCULATION /\/\" + print *, "Source batches: ", numToChar(self % N_steps) + print *, "Population per batch: ", numToChar(self % pop) + print *, "Initial RNG Seed: ", numToChar(self % pRNG % getSeed()) + print * + print *, repeat("<>",50) + end subroutine printSettings + +end module IMCPhysicsPackage_class diff --git a/PhysicsPackages/physicsPackageFactory_func.f90 b/PhysicsPackages/physicsPackageFactory_func.f90 index 9cda76b06..5d82bddbb 100644 --- a/PhysicsPackages/physicsPackageFactory_func.f90 +++ b/PhysicsPackages/physicsPackageFactory_func.f90 @@ -15,7 +15,7 @@ module physicsPackageFactory_func use fixedSourcePhysicsPackage_class, only : fixedSourcePhysicsPackage use vizPhysicsPackage_class, only : vizPhysicsPackage use rayVolPhysicsPackage_class, only : rayVolPhysicsPackage -! use dynamPhysicsPackage_class, only : dynamPhysicsPackage + use IMCPhysicsPackage_class, only : IMCPhysicsPackage implicit none private @@ -27,6 +27,7 @@ module physicsPackageFactory_func ! For now it is necessary to adjust trailing blanks so all enteries have the same length character(nameLen),dimension(*),parameter :: AVAILABLE_physicsPackages = [ 'eigenPhysicsPackage ',& 'fixedSourcePhysicsPackage',& + 'IMCPhysicsPackage ',& 'vizPhysicsPackage ',& 'rayVolPhysicsPackage '] @@ -70,14 +71,14 @@ function new_physicsPackage(dict) result(new) type is (fixedSourcePhysicsPackage) call new % init(dict) end select -! -! case('dynamPhysicsPackage') -! ! Allocate and initialise -! allocate( dynamPhysicsPackage :: new) -! select type(new) -! type is (dynamPhysicsPackage) -! call new % init(dict) -! end select + + case('IMCPhysicsPackage') + ! Allocate and initialise + allocate( IMCPhysicsPackage :: new) + select type(new) + type is (IMCPhysicsPackage) + call new % init(dict) + end select case('vizPhysicsPackage') diff --git a/SharedModules/CMakeLists.txt b/SharedModules/CMakeLists.txt index 8e6be469b..fd028b8a4 100644 --- a/SharedModules/CMakeLists.txt +++ b/SharedModules/CMakeLists.txt @@ -12,6 +12,7 @@ add_sources( ./genericProcedures.f90 ./statisticalTests_func.f90 ./timer_mod.f90 ./charLib_func.f90 + ./poly_func.f90 ./openmp_func.f90) add_unit_tests( ./Tests/grid_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/poly_func_test.f90) diff --git a/SharedModules/Tests/poly_func_test.f90 b/SharedModules/Tests/poly_func_test.f90 new file mode 100644 index 000000000..ef845e70f --- /dev/null +++ b/SharedModules/Tests/poly_func_test.f90 @@ -0,0 +1,46 @@ +module poly_func_test + use numPrecision + use poly_func, only : poly_integrate, poly_eval, poly_solve + use pFUnit_mod + + implicit none + +contains + +@Test + subroutine testPoly() + real(defReal), dimension(6) :: poly1, poly2, poly3, poly4 + real(defReal) :: x1, x2 + real(defReal) :: tol = 1.0E-4_defReal + + ! Test array + poly1(:) = [23.20_defReal, 1.59_defReal, 0.12_defReal, & + 10.02_defReal, 0.06_defReal, 8.03_defReal ] + + ! Analytical integral + poly2(:) = [2.10526_defReal, 1.5_defReal, 0.013289_defReal, & + 11.02_defReal, 1.06_defReal, 9.03_defReal ] + + ! Integrate using poly_integrate + poly3 = poly1 + call poly_integrate(poly3) + + ! Check + @assertEqual(poly3,poly2,tol) + + ! Evaluate using poly_eval + x1 = poly_eval(poly1, 1.074_defReal) + + ! Check + @assertEqual(x1,49.2504_defReal,tol) + + ! Solve poly2 = 1 using poly_solve, using x0 = 10 + x2 = poly_solve(poly2, poly1, 10.0_defReal, 1.0_defReal) + + ! Check + @assertEqual(x2,0.66643669_defReal,tol) + + end subroutine testPoly + + +end module poly_func_test diff --git a/SharedModules/endfConstants.f90 b/SharedModules/endfConstants.f90 index cbef7d1e6..e2e6d9b60 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -106,6 +106,7 @@ module endfConstants ! List of Macro MT numbers for macroscopic XSs. Unique to SCONE (not from Serpent) integer(shortInt), parameter :: macroAllScatter = -20 ,& macroAbsorbtion = -21 ,& + macroPlanck = -22 ,& noInteraction = -901 diff --git a/SharedModules/poly_func.f90 b/SharedModules/poly_func.f90 new file mode 100644 index 000000000..166f32f69 --- /dev/null +++ b/SharedModules/poly_func.f90 @@ -0,0 +1,165 @@ +module poly_func + + use universalVariables + use numPrecision + use genericProcedures + + implicit none + + !! Module to perform operations on simple polynomials + !! Polynomials are given as a 1D array containing n coefficients followed by n exponents + !! + !! Interface: + !! poly_integrate -> Update exponents and coefficients to integrate polynomial + !! - Currently gives error for input exponent of -1 + !! poly_solve -> Solves polynomial using Newton-Raphson method + !! + !! poly_eval -> Evaluates polynomial at given value + !! + + contains + + !! + !! Integrates a simple polynomial given coefficients and exponents, + !! returning indefinite integral + !! + !! Args: + !! equation -> 1D array of n coefficients followed by n exponents + !! + !! Errors: + !! Input array size is not divisible by 2 + !! Exponent of -1 is integrated + !! + subroutine poly_integrate(equation) + real(defReal), dimension(:), intent(inout) :: equation + integer(shortInt) :: n, i + character(100), parameter :: Here = "poly_integrate (poly_func.f90)" + + ! Check that array is of even size + if( modulo(size(equation), 2) /= 0 ) then + call fatalError(Here, "Array size must be divisible by 2") + end if + + n = size(equation) / 2 + + ! Integrate + do i=1, n + ! Update exponents + equation(i+n) = equation(i+n) + 1 + if( equation(i+n) == 0 ) then + call fatalError(Here, "Integrating exponent of -1, currently not yet supported") + end if + ! Update coefficients + equation(i) = equation(i) / equation(i+n) + end do + + end subroutine poly_integrate + + + !! + !! Use Newton-Raphspon method to solve polynomial with m terms + !! + !! Args: + !! equation -> 1D array of n coefficients followed by m exponents + !! derivative -> 1D array of n coefficients followed by m exponents + !! x0 -> Starting guess + !! const -> For f(x) = const, if not given then solves f(x) = 0 + !! + !! Errors: + !! Equation and derivative are different sizes + !! Input array sizes are not divisible by 2 + !! + function poly_solve(equation, derivative, x0, const) result(x) + real(defReal), dimension(:), intent(in) :: equation + real(defReal), dimension(:), intent(in) :: derivative + real(defReal), intent(in) :: x0 + real(defReal), intent(in), optional :: const + real(defReal) :: x, x_old, f, f_dash, c, tol + integer(shortInt) :: i, j, m + character(100), parameter :: Here = "poly_solve (poly_func.f90)" + + ! Check that input arrays are of same size + if( size(equation) /= size(derivative) ) then + call fatalError(Here, "Equation and Derivative array inputs are of different sizes") + end if + + ! Check that array is of even size + if( modulo(size(equation), 2) /= 0 ) then + call fatalError(Here, "Array size must be divisible by 2") + end if + + x = x0 + m = size(equation) / 2 + + ! May not converge if x0 = 0 + if ( x == 0 ) x = 0.0000001 + + ! If no constant present then solving f(x) = 0 + if( present(const) ) then + c = const + else + c = 0 + end if + + ! Iterate + i = 0 + iterate: do + ! Store x for convergence check + x_old = x + + ! Calculate f(x) and f'(x) + f = -c + f_dash = 0 + do j=1,m + f = f + equation(j) * x ** equation(j+m) + f_dash = f_dash + derivative(j) * x ** derivative(j+m) + end do + + ! Update x + x = x - f / f_dash + + ! Check for convergence + tol = 0.000000000001 + if( abs(x) > (1-tol)*abs(x_old) .and. abs(x) < (1+tol)*abs(x_old) ) exit iterate + + ! Call error if not converged + if( i >= 10000 ) then + call fatalError(Here, "Solution has not converged after 1000 iterations") + end if + + ! Increase counter + i = i+1 + + end do iterate + + end function poly_solve + + !! + !! Gives output value y for y = f(x) + !! + !! Args: + !! f -> Array defining polynomial + !! x -> Point at which to evaluate + !! + function poly_eval(f, x) result(y) + real(defReal), dimension(:), intent(in) :: f + real(defReal), intent(in) :: x + real(defReal) :: y + integer(shortInt) :: n, i + character(100), parameter :: Here = "poly_eval (poly_func.f90)" + + ! Check that array is of even size + if( modulo(size(f), 2) /= 0 ) then + call fatalError(Here, "Array size must be divisible by 2") + end if + + n = size(f) / 2 + + y = 0 + do i=1,n + y = y + f(i) * x ** f(i+n) + end do + + end function poly_eval + +end module poly_func diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index c62779f95..9f20f7a01 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -61,7 +61,9 @@ module universalVariables ! Particle Type Enumeration integer(shortInt), parameter :: P_NEUTRON_CE = 1, & - P_NEUTRON_MG = 2 + P_NEUTRON_MG = 2, & + P_PHOTON_CE = 3, & + P_PHOTON_MG = 4 ! Search error codes integer(shortInt), parameter :: valueOutsideArray = -1,& @@ -72,7 +74,10 @@ module universalVariables ! Physical constants real(defReal), parameter :: neutronMass = 939.5654133_defReal, & ! Neutron mass in MeV/c^2 lightSpeed = 2.99792458e10_defReal, & ! Light speed in cm/s - energyPerFission = 200.0_defReal ! MeV + energyPerFission = 200.0_defReal, & ! MeV + radiationConstant = 0.01372_defReal, & ! GJ/(cm^3 keV^4) + planckConst = 6.62607015e-30_defReal, & ! cm^2 kg/s + boltzmannConst = 1.380649e-19_defReal ! cm^2 kg s^-2 K^-1 ! Unit conversion real(defReal), parameter :: joulesPerMeV = 1.60218e-13 ! Convert MeV to J diff --git a/Tallies/TallyClerks/CMakeLists.txt b/Tallies/TallyClerks/CMakeLists.txt index f2d9ef595..6da7f659d 100644 --- a/Tallies/TallyClerks/CMakeLists.txt +++ b/Tallies/TallyClerks/CMakeLists.txt @@ -9,15 +9,16 @@ add_sources(./tallyClerk_inter.f90 ./keffImplicitClerk_class.f90 ./simpleFMClerk_class.f90 ./dancoffBellClerk_class.f90 - ./shannonEntropyClerk_class.f90 - ./centreOfMassClerk_class.f90 + ./shannonEntropyClerk_class.f90 + ./centreOfMassClerk_class.f90 + ./absorptionClerk_class.f90 ) add_unit_tests(./Tests/collisionClerk_test.f90 ./Tests/trackClerk_test.f90 ./Tests/keffAnalogClerk_test.f90 ./Tests/keffImplicitClerk_test.f90 - ./Tests/shannonEntropyClerk_test.f90 + ./Tests/shannonEntropyClerk_test.f90 ./Tests/simpleFMClerk_test.f90 ./Tests/collisionProbabilityClerk_test.f90 ) diff --git a/Tallies/TallyClerks/absorptionClerk_class.f90 b/Tallies/TallyClerks/absorptionClerk_class.f90 new file mode 100644 index 000000000..9426b63bf --- /dev/null +++ b/Tallies/TallyClerks/absorptionClerk_class.f90 @@ -0,0 +1,322 @@ +module absorptionClerk_class + + use numPrecision + use tallyCodes + use genericProcedures, only : fatalError + use dictionary_class, only : dictionary + use particle_class, only : particle, particleState + use outputFile_class, only : outputFile + use scoreMemory_class, only : scoreMemory + use tallyClerk_inter, only : tallyClerk, kill_super => kill + + ! Nuclear Data interface + use nuclearDatabase_inter, only : nuclearDatabase + + ! Tally Filters + use tallyFilter_inter, only : tallyFilter + use tallyFilterFactory_func, only : new_tallyFilter + + ! Tally Maps + use tallyMap_inter, only : tallyMap + use tallyMapFactory_func, only : new_tallyMap + + ! Tally Responses + use tallyResponseSlot_class, only : tallyResponseSlot + + + use tallyResult_class, only : tallyResult, tallyResultEmpty + + implicit none + private + + !! + !! Colision estimator of reaction rates + !! Calculates flux weighted integral from collisions + !! + !! Private Members: + !! filter -> Space to store tally Filter + !! map -> Space to store tally Map + !! response -> Array of responses + !! width -> Number of responses (# of result bins for each map position) + !! + !! Interface + !! tallyClerk Interface + !! + !! SAMPLE DICTIOANRY INPUT: + !! + !! myAbsorptionClerk { + !! type absorptionClerk; + !! # filter { } # + !! # map { } # + !! response (resName1 #resName2 ... #) + !! resName1 { } + !! #resNamew { run-time procedures + procedure :: reportHist + + ! Output procedures + procedure :: display + procedure :: print + procedure :: getResult + + end type absorptionClerk + + + type,public, extends(tallyResult) :: absClerkResult + real(defReal), dimension(:), allocatable :: clerkResults + end type absClerkResult + + +contains + + !! + !! Initialise clerk from dictionary and name + !! + !! See tallyClerk_inter for details + !! + subroutine init(self, dict, name) + class(absorptionClerk), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(nameLen), intent(in) :: name + character(nameLen),dimension(:),allocatable :: responseNames + integer(shortInt) :: i + + ! Assign name + call self % setName(name) + + ! Load filter + if( dict % isPresent('filter')) then + call new_tallyFilter(self % filter, dict % getDictPtr('filter')) + end if + + ! Load map + if( dict % isPresent('map')) then + call new_tallyMap(self % map, dict % getDictPtr('map')) + end if + + ! Get names of response dictionaries + call dict % get(responseNames,'response') + + ! Load responses + allocate(self % response(size(responseNames))) + do i=1, size(responseNames) + call self % response(i) % init(dict % getDictPtr( responseNames(i) )) + end do + + ! Set width + self % width = size(responseNames) + + end subroutine init + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(absorptionClerk), intent(inout) :: self + + ! Superclass + call kill_super(self) + + ! Kill and deallocate filter + if(allocated(self % filter)) then + deallocate(self % filter) + end if + + ! Kill and deallocate map + if(allocated(self % map)) then + call self % map % kill() + deallocate(self % map) + end if + + ! Kill and deallocate responses + if(allocated(self % response)) then + deallocate(self % response) + end if + + self % width = 0 + + end subroutine kill + + !! + !! Returns array of codes that represent diffrent reports + !! + !! See tallyClerk_inter for details + !! + function validReports(self) result(validCodes) + class(absorptionClerk),intent(in) :: self + integer(shortInt),dimension(:),allocatable :: validCodes + + validCodes = [hist_CODE] + + end function validReports + + !! + !! Return memory size of the clerk + !! + !! See tallyClerk_inter for details + !! + elemental function getSize(self) result(S) + class(absorptionClerk), intent(in) :: self + integer(shortInt) :: S + + S = size(self % response) + if(allocated(self % map)) S = S * self % map % bins(0) + + end function getSize + + !! + !! Process incoming collision report + !! + !! See tallyClerk_inter for details + !! + subroutine reportHist(self, p, xsData, mem) + class(absorptionClerk), intent(inout) :: self + class(particle), intent(in) :: p + class(nuclearDatabase), intent(inout) :: xsData + type(scoreMemory), intent(inout) :: mem + type(particleState) :: state + integer(shortInt) :: binIdx, i + integer(longInt) :: adrr + real(defReal) :: scoreVal, flx + character(100), parameter :: Here =' reportHist (absorptionClerk_class.f90)' + + ! Get current particle state + state = p + + if( p % fate == LEAK_FATE ) return + if( p % isDead .eqv. .false.) call fatalError(Here, 'Particle is still alive') + + ! Check if within filter + if(allocated( self % filter)) then + if(self % filter % isFail(state)) return + end if + + ! Find bin index + if(allocated(self % map)) then + binIdx = self % map % map(state) + else + binIdx = 1 + end if + + ! Return if invalid bin index + if (binIdx == 0) return + + ! Calculate bin address + adrr = self % getMemAddress() + self % width * (binIdx -1) - 1 + + ! Calculate flux sample 1/totXs + flx = ONE / xsData % getTotalMatXS(p, p % matIdx()) + + ! Append all bins + do i=1,self % width + scoreVal = self % response(i) % get(p, xsData) * p % w *flx + call mem % score(scoreVal, adrr + i) + + end do + + end subroutine reportHist + + !! + !! Display convergance progress on the console + !! + !! See tallyClerk_inter for details + !! + subroutine display(self, mem) + class(absorptionClerk), intent(in) :: self + type(scoreMemory), intent(in) :: mem + + print *, 'absorptionClerk does not support display yet' + + end subroutine display + + !! + !! Write contents of the clerk to output file + !! + !! See tallyClerk_inter for details + !! + subroutine print(self, outFile, mem) + class(absorptionClerk), intent(in) :: self + class(outputFile), intent(inout) :: outFile + type(scoreMemory), intent(in) :: mem + real(defReal) :: val, std + integer(shortInt) :: i + integer(shortInt),dimension(:),allocatable :: resArrayShape + character(nameLen) :: name + + ! Begin block + call outFile % startBlock(self % getName()) + + ! If collision clerk has map print map information + if( allocated(self % map)) then + call self % map % print(outFile) + end if + + ! Write results. + ! Get shape of result array + if(allocated(self % map)) then + resArrayShape = [size(self % response), self % map % binArrayShape()] + else + resArrayShape = [size(self % response)] + end if + + ! Start array + name ='Res' + call outFile % startArray(name, resArrayShape) + + ! Print results to the file + do i=1,product(resArrayShape) + call mem % getResult(val, std, self % getMemAddress() - 1 + i) + call outFile % addResult(val,std) + + end do + + call outFile % endArray() + call outFile % endBlock() + + end subroutine print + + !! + !! Return result for interaction with Physics Package + !! from the clerk in the slot + !! + !! See tallyClerk_inter for details + !! + pure subroutine getResult(self, res, mem) + class(absorptionClerk), intent(in) :: self + class(tallyResult), allocatable, intent(inout) :: res + type(scoreMemory), intent(in) :: mem + real(defReal), dimension(:), allocatable :: w + integer(shortInt) :: i, N + + N = self % getSize() + allocate( w(N) ) + + ! Get result value for each material + do i = 1, N + call mem % getResult(w(i), self % getMemAddress()+i-1) + end do + + allocate(res, source = absClerkResult(w)) + + end subroutine getResult + +end module absorptionClerk_class diff --git a/Tallies/TallyClerks/collisionClerk_class.f90 b/Tallies/TallyClerks/collisionClerk_class.f90 index 242fceb39..ef7f8a7a1 100644 --- a/Tallies/TallyClerks/collisionClerk_class.f90 +++ b/Tallies/TallyClerks/collisionClerk_class.f90 @@ -93,7 +93,7 @@ subroutine init(self, dict, name) ! Assign name call self % setName(name) - ! Load filetr + ! Load filter if( dict % isPresent('filter')) then call new_tallyFilter(self % filter, dict % getDictPtr('filter')) end if diff --git a/Tallies/TallyClerks/tallyClerkFactory_func.f90 b/Tallies/TallyClerks/tallyClerkFactory_func.f90 index b3bbced14..bc4579273 100644 --- a/Tallies/TallyClerks/tallyClerkFactory_func.f90 +++ b/Tallies/TallyClerks/tallyClerkFactory_func.f90 @@ -17,6 +17,7 @@ module tallyClerkFactory_func use dancoffBellClerk_class, only : dancoffBellClerk use shannonEntropyClerk_class, only : shannonEntropyClerk use centreOfMassClerk_class, only : centreOfMassClerk + use absorptionClerk_class, only : absorptionClerk implicit none private @@ -36,7 +37,8 @@ module tallyClerkFactory_func 'simpleFMClerk ',& 'shannonEntropyClerk ',& 'centreOfMassClerk ',& - 'dancoffBellClerk '] + 'dancoffBellClerk ',& + 'absorptionClerk ' ] contains @@ -96,6 +98,10 @@ subroutine new_tallyClerk(new, dict, name) allocate(centreOfMassClerk :: new) call new % init(dict, name) + case('absorptionClerk') + allocate(absorptionClerk :: new) + call new % init(dict, name) + !*** NEW TALLY MAP TEMPLATE ***! !case('') ! allocate( :: new) diff --git a/Tallies/TallyClerks/tallyClerk_inter.f90 b/Tallies/TallyClerks/tallyClerk_inter.f90 index e73904f66..205cc5cd9 100644 --- a/Tallies/TallyClerks/tallyClerk_inter.f90 +++ b/Tallies/TallyClerks/tallyClerk_inter.f90 @@ -47,7 +47,7 @@ module tallyClerk_inter !! Interface: !! init -> Initialise Clerk from a dictionary !! kill -> Return to Uninitialised State - !! validReports -> Returns array of integers with tallyCodes of reports that the clark requires + !! validReports -> Returns array of integers with tallyCodes of reports that the clerk requires !! getSize -> Return size required by Clerk on ScoreMemory !! setMemAddress -> Setter for "memAddress" member !! getMemAddress -> Getter for "memAddress" member diff --git a/Tallies/TallyResponses/tallyResponse_inter.f90 b/Tallies/TallyResponses/tallyResponse_inter.f90 index 314d6c0e0..d677b100d 100644 --- a/Tallies/TallyResponses/tallyResponse_inter.f90 +++ b/Tallies/TallyResponses/tallyResponse_inter.f90 @@ -21,10 +21,10 @@ module tallyResponse_inter !! !! Interface: !! init -> Initialise - !! get -> Get velue of the response + !! get -> Get value of the response !! kill -> Return to uninitialised state !! - type, public,abstract :: tallyResponse + type, public, abstract :: tallyResponse private contains diff --git a/Tallies/TallyResponses/weightResponse_class.f90 b/Tallies/TallyResponses/weightResponse_class.f90 index e007e20dd..87c599564 100644 --- a/Tallies/TallyResponses/weightResponse_class.f90 +++ b/Tallies/TallyResponses/weightResponse_class.f90 @@ -4,12 +4,14 @@ module weightResponse_class use endfConstants use genericProcedures, only : fatalError, numToChar use dictionary_class, only : dictionary - use particle_class, only : particle, P_NEUTRON + use particle_class, only : particle use tallyResponse_inter, only : tallyResponse ! Nuclear Data interfaces use nuclearDatabase_inter, only : nuclearDatabase - use neutronMaterial_inter, only : neutronMaterial, neutronMaterial_CptrCast + use materialHandle_inter, only : materialHandle + use neutronMaterial_inter, only : neutronMaterial_CptrCast + use imcMaterial_inter, only : IMCMaterial_CptrCast implicit none private @@ -70,17 +72,15 @@ function get(self, p, xsData) result(val) class(particle), intent(in) :: p class(nuclearDatabase), intent(inout) :: xsData real(defReal) :: val - class(neutronMaterial), pointer :: mat + class(materialHandle), pointer :: mat val = ZERO - ! Return 0.0 if particle is not neutron - if(p % type /= P_NEUTRON) return - ! Get pointer to active material data mat => neutronMaterial_CptrCast(xsData % getMaterial(p % matIdx())) + if(.not.associated(mat)) mat => IMCMaterial_CptrCast(xsData % getMaterial(p % matIdx())) - ! Return if material is not a neutronMaterial + ! Return if material is not a neutronMaterial or IMCMaterial if(.not.associated(mat)) return if (self % moment == 0) then diff --git a/Tallies/scoreMemory_class.f90 b/Tallies/scoreMemory_class.f90 index 6b8d15a54..a309c6965 100644 --- a/Tallies/scoreMemory_class.f90 +++ b/Tallies/scoreMemory_class.f90 @@ -53,12 +53,14 @@ module scoreMemory_class !! cumulative sums. Then sets the bin to zero. !! !! closeCycle(normFactor): Multiplies all scores by normFactor and accumulates them in - !! cumulative sums. Sets all scors to zero. + !! cumulative sums. Sets all scores to zero. !! !! lastCycle(): Return true if the next call to closeCycle will close a batch. !! !! getBatchSize(): Returns number of cycles that constitute a single batch. !! + !! reset(idx): Reset the value in a memory location + !! !! Example use case: !! !! do batches=1,20 @@ -99,6 +101,7 @@ module scoreMemory_class procedure :: closeBin procedure :: lastCycle procedure :: getBatchSize + procedure :: reset ! Private procedures procedure, private :: score_defReal @@ -276,7 +279,6 @@ subroutine closeCycle(self, normFactor) ! Increment Cycle Counter self % cycles = self % cycles + 1 - if(mod(self % cycles, self % batchSize) == 0) then ! Close Batch !$omp parallel do @@ -446,4 +448,19 @@ elemental function getScore(self, idx) result (score) end function getScore + !! + !! Reset the value in a memory location + !! Useful for tallies that do not accumulate across cycles + !! + subroutine reset(self,idx) + class(scoreMemory), intent(inout) :: self + integer(longInt), intent(in) :: idx + + self % bins(idx, :) = ZERO + + self % cycles = 0 + self % batchN = 0 + + end subroutine reset + end module scoreMemory_class diff --git a/Tallies/tallyAdmin_class.f90 b/Tallies/tallyAdmin_class.f90 index 067a52c3d..f5cf813e6 100644 --- a/Tallies/tallyAdmin_class.f90 +++ b/Tallies/tallyAdmin_class.f90 @@ -73,6 +73,7 @@ module tallyAdmin_class !! 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 + !! reset -> Resets tally clerk count to 0 !! !! SAMPLE DICTIOANRY INPUT: !! @@ -137,6 +138,7 @@ module tallyAdmin_class ! Interaction procedures procedure :: getResult + procedure :: reset ! Display procedures procedure :: display @@ -747,6 +749,37 @@ pure subroutine getResult(self, res, name) end subroutine getResult + !! + !! Resets tally clerk count to 0 + !! + subroutine reset(self, name) + class(tallyAdmin),intent(inout) :: self + character(*), intent(in) :: name + character(nameLen) :: name_loc + integer(shortInt) :: idx + integer(shortInt),parameter :: NOT_PRESENT = -3 + integer(longInt) :: addr + integer(shortInt) :: i, width + character(100),parameter :: Here='reset (tallyAdmin_class.f90)' + + name_loc = name + + idx = self % clerksNameMap % getOrDefault(name_loc, NOT_PRESENT) + + if(idx == NOT_PRESENT) then + call fatalError(Here, 'Tally clerk not present') + end if + + addr = self % tallyClerks(idx) % getMemAddress() + width = self % tallyClerks(idx) % getSize() + + do i = 1, width + call self % mem % reset(addr+i-1) + end do + + end subroutine reset + + !! !! Append sorting array identified with the code with tallyClerk idx !! diff --git a/Tallies/tallyCodes.f90 b/Tallies/tallyCodes.f90 index 44f7bac69..518696c20 100644 --- a/Tallies/tallyCodes.f90 +++ b/Tallies/tallyCodes.f90 @@ -22,6 +22,7 @@ module tallyCodes integer(shortInt),parameter,public :: abs_FATE = 5000 ,& leak_FATE = 5001 ,& lost_FATE = 5002 ,& - aged_FATE = 5003 + aged_FATE = 5003 ,& + time_FATE = 5004 end module tallyCodes diff --git a/TransportOperator/CMakeLists.txt b/TransportOperator/CMakeLists.txt index 2d71564c0..259e709e5 100644 --- a/TransportOperator/CMakeLists.txt +++ b/TransportOperator/CMakeLists.txt @@ -4,4 +4,6 @@ add_sources(./transportOperator_inter.f90 ./transportOperatorDT_class.f90 # ./transportOperatorDynamicDT_class.f90 ./transportOperatorST_class.f90 - ./transportOperatorHT_class.f90) + ./transportOperatorHT_class.f90 + ./transportOperatorTimeHT_class.f90 + ./Grid/trackingGrid_class.f90) diff --git a/TransportOperator/Grid/trackingGrid_class.f90 b/TransportOperator/Grid/trackingGrid_class.f90 new file mode 100644 index 000000000..aa9a4a4de --- /dev/null +++ b/TransportOperator/Grid/trackingGrid_class.f90 @@ -0,0 +1,403 @@ +module trackingGrid_class + + use numPrecision + use universalVariables, only : SURF_TOL, P_PHOTON_MG + use genericProcedures, only : fatalError, numToChar + use dictionary_class, only : dictionary + use geometry_inter, only : geometry + use dynArray_class, only : dynIntArray + use nuclearDatabase_inter, only : nuclearDatabase + use particle_class, only : particle + + use geometryReg_mod, only : gr_geomPtr => geomPtr + use nuclearDataReg_mod, only : ndReg_get => get + + !! + !! + !! + type, private :: gridCell + integer(shortInt), dimension(:), allocatable :: mats + real(defReal) :: majorant + + end type gridCell + + !! + !! As in latUniverse_class, idx is 1 in bottom X, Y & Z corner. + !! It increases first with X then Y and lastly Z. + !! + !! sizeN -> array [nx, ny, nz], the dimensions of the grid + !! pitch -> array [dx, dy, dz], the discretisation in each direction + !! bounds -> [x_min, y_min, z_min, z_max, y_max, z_max] as in geometry_inter + !! + type, public :: trackingGrid + class(geometry), pointer :: mainGeom => null() + class(nuclearDatabase), pointer :: xsData => null() + integer(shortInt), dimension(:), allocatable :: sizeN + real(defReal), dimension(3) :: pitch = 0 + real(defReal), dimension(6) :: bounds + real(defReal), dimension(3) :: corner + real(defReal), dimension(3) :: a_bar + type(gridCell), dimension(:), allocatable :: gridCells + + contains + procedure :: init + !procedure :: kill + procedure :: getDistance + procedure :: getValue + procedure :: storeMats + procedure :: update + + end type trackingGrid + +contains + + subroutine init(self, dict, geom, xsData) + class(trackingGrid), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), intent(in), pointer, optional :: geom + class(nuclearDatabase), intent(in), pointer, optional :: xsData + integer(shortInt) :: N + integer(shortInt), dimension(:), allocatable :: searchN + + ! Store pointer to main geometry and data +! self % mainGeom => geom +! self % xsData => xsData + self % xsData => ndReg_get(P_PHOTON_MG) ! TODO: not an ideal way to do this but fine temporarily + self % mainGeom => gr_geomPtr(1) + + ! Store settings + call dict % get(self % sizeN, 'dimensions') + call dict % get(searchN, 'searchN') + + ! Get bounds of grid and calculate discretisations + self % bounds = self % mainGeom % bounds() + + self % pitch(1) = (self % bounds(4) - self % bounds(1)) / self % sizeN(1) + self % pitch(2) = (self % bounds(5) - self % bounds(2)) / self % sizeN(2) + self % pitch(3) = (self % bounds(6) - self % bounds(3)) / self % sizeN(3) + + self % corner = [self % bounds(1), self % bounds(2), self % bounds(3)] + self % a_bar = self % pitch * HALF - SURF_TOL + + ! Allocate space for cells + N = self % sizeN(1) * self % sizeN(2) * self % sizeN(3) + allocate(self % gridCells(N)) + + ! Find material idxs present in each cell + call self % storeMats(searchN) + + end subroutine init + + + !! + !! May have issues with non-box geometry root universe surface with reflective boundary + !! + function getDistance(self, r, u) result(dist) + class(trackingGrid), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal) :: dist + real(defReal), dimension(3) :: r_bar, low, high !, point, corner, ratio + character(100), parameter :: Here = 'getDistance (trackingGrid_class.f90)' + + ! Calculate position from grid corner + r_bar = r - self % corner + if (any(r_bar < -SURF_TOL)) call fatalError(Here, 'Point is outside grid geometry') !TODO only checks bottom for now + + ! Write as a fraction across cell + r_bar = r_bar / self % pitch + r_bar = r_bar - floor(r_bar) + + ! Account for surface tolerance + low = SURF_TOL / self % pitch + high = ONE - low + do i = 1, 3 + if (r_bar(i) < low(i) .and. u(i) < ZERO) r_bar(i) = ONE + if (r_bar(i) > high(i) .and. u(i) > ZERO) r_bar(i) = ZERO + end do + + ! Distance to centre plus distance from centre to required boundary + r_bar = (HALF - r_bar + sign(HALF, u)) * self % pitch + dist = minval(r_bar / u) + + if (dist <= ZERO) call fatalError(Here, 'Distance invalid: '//numToChar(dist)) + + ! Increase by surface tolerance to ensure that boundary conditions are correctly applied + dist = dist + SURF_TOL + + + ! Round each dimension either up or down depending on which boundary will be hit +! do i = 1, 3 +! if (u(i) >= 0) then +! ! Round each dimension either up or down depending on which boundary will be hit +! do i = 1, 3 +! if (u(i) >= 0) then +! corner(i) = ceiling(point(i)) +! else +! corner(i) = floor(point(i)) +! end if +! ! Adjust if starting position was on boundary +! if (abs(corner(i) - point(i)) < SURF_TOL) then +! corner(i) = corner(i) + sign(ONE, u(i)) +! end if +! end do + +! ! Convert back to spatial coordinates - this is now the coordinates of the corner being travelled towards +! corner = corner * self % pitch + +! ! Determine which axis boundary will be hit first +! ratio = (corner - r_bar) / u + +! dist = minval(ratio) + + + end function getDistance + + + !! + !! Returns value of grid cell at position + !! + function getValue(self, r, u) result(val) + class(trackingGrid), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal) :: val + real(defReal), dimension(3) :: r_bar + integer(shortInt), dimension(3) :: corner, ijk + integer(shortInt) :: i, idx + character(100), parameter :: Here = 'getValue (trackingGrid_class.f90)' + + ! Find lattice location in x,y&z + ijk = floor((r - self % corner) / self % pitch) + 1 + + ! Get position wrt middle of the lattice cell + r_bar = r - self % corner - ijk * self % pitch + HALF * self % pitch + + ! Check if position is within surface tolerance + ! If it is, push it to next cell + do i = 1, 3 + if (abs(r_bar(i)) > self % a_bar(i) .and. r_bar(i)*u(i) > ZERO) then + + ! Select increment. Ternary expression + if (u(i) < ZERO) then + inc = -1 + else + inc = 1 + end if + + ijk(i) = ijk(i) + inc + + end if + end do + + ! Set localID & cellIdx + if (any(ijk <= 0 .or. ijk > self % sizeN)) then ! Point is outside grid + call fatalError(Here, 'Point is outside grid') + + else + idx = ijk(1) + self % sizeN(1) * (ijk(2)-1 + self % sizeN(2) * (ijk(3)-1)) + + end if + + + +! ! Get grid cell bottom corner +! r_bar = reposition(r, self % bounds) - self % corner +! corner = floor(r_bar) +! do i = 1, 3 +! if (corner(i) == r_bar(i) .and. u(i) < 0) then +! ! Adjust for point starting on cell boundary +! corner(i) = corner(i) - 1 +! end if +! end do +! +! ! Adjust for bottom corner starting at 1 +! corner = corner + 1 +! +! ! Get grid cell idx +! idx = get_idx(corner, self % sizeN) +! if (idx == 0) call fatalError(Here, 'Point is outside grid: '//numToChar(r)) + + val = self % gridCells(idx) % majorant + + if (val < ZERO) call fatalError(Here, 'Invalid majorant: '//numToChar(val)) + + end function getValue + + !! + !! + !! + subroutine storeMats(self, searchN) + class(trackingGrid), intent(inout) :: self + integer(shortInt), dimension(3), intent(in) :: searchN + real(defReal), dimension(3) :: searchRes + integer(shortInt) :: i, j, k, l, matIdx, uniqueID + real(defReal), dimension(3) :: corner, r + type(dynIntArray) :: mats + + ! Calculate distance between search points + searchRes = self % pitch / (searchN + 1) + + ! Loop through grid cells + do i = 1, size(self % gridCells) + + ! Get cell lower corner + corner = self % corner + self % pitch * (get_ijk(i, self % sizeN) - 1) + + ! Loop through search locations + do j = 1, searchN(1) + do k = 1, searchN(2) + do l = 1, searchN(3) + ! Find matIdx at search location + r = corner + [j, k, l] * searchRes + call self % mainGeom % whatIsAt(matIdx, uniqueID, r) + + ! Add to array if not already present + if (mats % isPresent(matIdx)) then + ! Do nothing + else + call mats % add(matIdx) + end if + + end do + end do + end do + + ! Store matIdx data in grid cell + self % gridCells(i) % mats = mats % expose() + call mats % kill() + + end do + + end subroutine storeMats + + !! + !! + !! + subroutine update(self) + class(trackingGrid), intent(inout) :: self + integer(shortInt) :: i + integer(shortInt), save :: j, matIdx + real(defReal), save :: sigmaT + class(particle), allocatable :: p + !$omp threadprivate(j, matIdx) + + allocate(p) + p % G = 1 + + !$omp parallel do + ! Loop through grid cells + do i = 1, size(self % gridCells) + ! Reset majorant + self % gridCells(i) % majorant = ZERO + + do j = 1, size(self % gridCells(i) % mats) + ! Get opacity of each material + matIdx = self % gridCells(i) % mats(j) + if (matIdx /= 0) then + sigmaT = self % xsData % getTransMatXS(p, matIdx) + ! Update majorant if required + if (sigmaT > self % gridCells(i) % majorant) self % gridCells(i) % majorant = sigmaT + end if + + end do + end do + !$omp end parallel do + + end subroutine update + + + + !! + !! Generate ijk from localID and shape + !! + !! Args: + !! localID [in] -> Local id of the cell between 1 and product(sizeN) + !! sizeN [in] -> Number of cells in each cardinal direction x, y & z + !! + !! Result: + !! Array ijk which has integer position in each cardinal direction + !! + pure function get_ijk(localID, sizeN) result(ijk) + integer(shortInt), intent(in) :: localID + integer(shortInt), dimension(3), intent(in) :: sizeN + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: temp, base + + temp = localID - 1 + + base = temp / sizeN(1) + ijk(1) = temp - sizeN(1) * base + 1 + + temp = base + base = temp / sizeN(2) + ijk(2) = temp - sizeN(2) * base + 1 + + ijk(3) = base + 1 + + end function get_ijk + + + pure function get_idx(ijk, sizeN) result(idx) + integer(shortInt), dimension(3), intent(in) :: ijk + integer(shortInt), dimension(3), intent(in) :: sizeN + integer(shortInt) :: idx + + if (any(ijk <= 0 .or. ijk > sizeN)) then ! Point is outside grid + idx = 0 + else + idx = ijk(1) + sizeN(1) * (ijk(2)-1 + sizeN(2) * (ijk(3)-1)) + end if + + end function get_idx + + !! + !! Adjustment for surface tolerance used by getValue subroutine + !! + function reposition(r, bounds) result(rNew) + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(6), intent(in) :: bounds + real(defReal), dimension(3) :: rNew + integer(shortInt) :: i + + rNew = r + + do i = 1, 3 + if (r(i) < bounds(i) .and. r(i) > bounds(i) -SURF_TOL) rNew(i) = bounds(i) + if (r(i) > bounds(i+3) .and. r(i) < bounds(i+3)+SURF_TOL) rNew(i) = bounds(i+3) + end do + + ! TODO Boundaries between cells rather than just edge of grid + + end function reposition + + !! + !! Adjustment for surface tolerance used by getDistance function. + !! Able to be simpler than repositionLoc as only consider position within cell + !! rather than within grid. + !! + !! Args: + !! r [inout] -> position as a fraction of distance across cell, 0 < r(i), < 1 + !! u [in] -> direction + !! pitch [in] -> grid resolution + !! +! subroutine repositionDist(r_bar, u, pitch) +! real(defReal), dimension(3), intent(inout) :: r_bar +! real(defReal), dimension(3), intent(in) :: u +! real(defReal), dimension(3), intent(in) :: pitch +! real(defReal), dimension(3) :: low, high +! integer(shortInt) :: i +! +! ! Calculate cut-offs +! low = SURF_TOL / pitch +! high = ONE - low +! +! ! Change position if needed +! do i = 1, 3 +! if (r_bar(i) < low(i) .and. u(i) < ZERO) r_bar(i) = ONE +! if (r_bar(i) > high(i) .and. u(i) > ZERO) r_bar(i) = ZERO +! end do +! +! end subroutine repositionDist + +end module trackingGrid_class diff --git a/TransportOperator/transportOperatorFactory_func.f90 b/TransportOperator/transportOperatorFactory_func.f90 index 67b821f95..2949c2956 100644 --- a/TransportOperator/transportOperatorFactory_func.f90 +++ b/TransportOperator/transportOperatorFactory_func.f90 @@ -12,6 +12,7 @@ module transportOperatorFactory_func use transportOperatorST_class, only : transportOperatorST use transportOperatorDT_class, only : transportOperatorDT use transportOperatorHT_class, only : transportOperatorHT + use transportOperatorTimeHT_class, only : transportOperatorTimeHT !use transportOperatorDynamicDT_class, only : transportOperatorDynamicDT implicit none @@ -22,9 +23,10 @@ module transportOperatorFactory_func ! It is printed if type was unrecognised ! NOTE: ! For now it is necessary to adjust trailing blanks so all enteries have the same length - character(nameLen),dimension(*),parameter :: AVALIBLE_transportOps = [ 'transportOperatorST', & - 'transportOperatorDT', & - 'transportOperatorHT']!, & + character(nameLen),dimension(*),parameter :: AVALIBLE_transportOps = [ 'transportOperatorST ', & + 'transportOperatorDT ', & + 'transportOperatorHT ', & + 'transportOperatorTimeHT']!, & ! 'dynamicTranspOperDT'] public :: new_transportOperator @@ -46,7 +48,7 @@ subroutine new_transportOperator(new, dict) ! Obtain string that specifies type to be built call dict % get(type,'type') - ! Allocate approperiate subclass of transportOperator + ! Allocate appropriate subclass of transportOperator ! *** ADD CASE STATEMENT FOR A NEW TRANSPORT OPERATOR BELOW ***! select case(type) case('transportOperatorST') @@ -61,6 +63,10 @@ subroutine new_transportOperator(new, dict) allocate( transportOperatorHT :: new) call new % init(dict) + case('transportOperatorTimeHT') + allocate( transportOperatorTimeHT :: new) + call new % init(dict) + ! case('dynamicTranspOperDT') ! allocate( transportOperatorDynamicDT :: new) ! call new % init(dict, geom) diff --git a/TransportOperator/transportOperatorTimeHT_class.f90 b/TransportOperator/transportOperatorTimeHT_class.f90 new file mode 100644 index 000000000..cd78f63d5 --- /dev/null +++ b/TransportOperator/transportOperatorTimeHT_class.f90 @@ -0,0 +1,284 @@ +!! +!! Transport operator for time-dependent problems using a hybrid of delta tracking and surface tracking +!! +module transportOperatorTimeHT_class + use numPrecision + use universalVariables + + use genericProcedures, only : fatalError, numToChar + use particle_class, only : particle, P_PHOTON + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary + use rng_class, only : rng + + ! Superclass + use transportOperator_inter, only : transportOperator, init_super => init + + ! Tally interface + use tallyCodes + use tallyAdmin_class, only : tallyAdmin + + ! Nuclear data interfaces + use nuclearDatabase_inter, only : nuclearDatabase + + implicit none + private + + !! + !! Tracking method + !! + integer(shortInt), parameter :: HT = 1 ! Hybrid tracking + integer(shortInt), parameter :: GT = 2 ! Grid tracking + integer(shortInt), parameter :: ST = 3 ! Surface tracking + integer(shortInt), parameter :: DT = 4 ! Delta tracking + + !! + !! Transport operator that moves a particle with using hybrid tracking, up to a time boundary + !! + type, public, extends(transportOperator) :: transportOperatorTimeHT + real(defReal) :: deltaT + real(defReal) :: cutoff + integer(shortInt) :: method + contains + procedure :: transit => timeTracking + procedure :: init + procedure, private :: surfaceTracking + procedure, private :: deltaTracking + procedure, private :: getMajInv + end type transportOperatorTimeHT + +contains + + subroutine timeTracking(self, p, tally, thisCycle, nextCycle) + class(transportOperatorTimeHT), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + class(particleDungeon), intent(inout) :: thisCycle + class(particleDungeon), intent(inout) :: nextCycle + real(defReal) :: sigmaT + character(100), parameter :: Here = 'timeTracking (transportOperatorTimeHT_class.f90)' + + ! Select action based on specified method - HT and GT start with DT but can switch to ST + if (self % method == ST) then + call self % surfaceTracking(p) + else + call self % deltaTracking(p) + end if + + ! Check for particle leakage + if (p % matIdx() == OUTSIDE_FILL) then + p % fate = LEAK_FATE + p % isDead = .true. + end if + + call tally % reportTrans(p) + + end subroutine timeTracking + + !! + !! Perform surface tracking + !! + subroutine surfaceTracking(self, p) + class(transportOperatorTimeHT), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime, dColl, dist, sigmaT + integer(shortInt) :: event + character(100), parameter :: Here = 'surfaceTracking (transportOperatorTimeHT_class.f90)' + + STLoop:do + + ! Find distance to time boundary + dTime = lightSpeed * (p % timeMax - p % time) + + ! Sample distance to collision + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + dColl = -log( p % pRNG % get() ) / sigmaT + + ! Choose minimum distance + dist = min(dTime, dColl) + + ! Move through geometry using minimum distance + call self % geom % move(p % coords, dist, event) + + ! Check for particle leakage + if (p % matIdx() == OUTSIDE_FILL) return + + ! Increase time based on distance moved + p % time = p % time + dist / lightSpeed + + ! Check result of transport + if (dist == dTime) then + ! Time boundary + if (event /= COLL_EV) call fatalError(Here, 'Move outcome should be COLL_EV & + &after moving dTime') + p % fate = AGED_FATE + p % time = p % timeMax + exit STLoop + + else if (dist == dColl) then + ! Collision, increase time accordingly + if (event /= COLL_EV) call fatalError(Here, 'Move outcome should be COLL_EV & + &after moving dTime') + exit STLoop + + end if + + end do STLoop + + end subroutine surfaceTracking + + !! + !! Perform delta tracking - option to switch to surface tracking for HT and GT methods + !! + subroutine deltaTracking(self, p) + class(transportOperatorTimeHT), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime, dColl, dGrid, sigmaT, majorant_inv, dist + character(100), parameter :: Here = 'deltaTracking (transportOperatorTimeHT_class.f90)' + + ! Get majorant and grid crossing distance if required + majorant_inv = self % getMajInv(p) + + ! Get initial opacity + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + + ! Check if surface tracking is needed, avoiding unnecessary grid calculations + if (sigmaT * majorant_inv < ONE - self % cutoff) then + call self % surfaceTracking(p) + return + end if + + ! Calculate initial distance to grid + if (self % method == GT) then + dGrid = self % grid % getDistance(p % coords % lvl(1) % r, p % coords % lvl(1) % dir) + else + dGrid = INF + end if + + DTLoop:do + + ! Find distance to time boundary + dTime = lightSpeed * (p % timeMax - p % time) + + ! Sample distance to collision + dColl = -log( p % pRNG % get() ) * majorant_inv + + ! Select particle by minimum distance + dist = min(dColl, dTime, dGrid) + call self % geom % teleport(p % coords, dist) + p % time = p % time + dist / lightSpeed + + ! Check for particle leakage + if (p % matIdx() == OUTSIDE_FILL) return + + ! Act based on distance moved + if (dist == dGrid) then + ! Update values and cycle loop + majorant_inv = self % getMajInv(p) + dGrid = self % grid % getDistance(p % coords % lvl(1) % r, p % coords % lvl(1) % dir) + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + cycle DTLoop + + else if (dist == dTime) then + ! Update particle fate and exit + p % fate = AGED_FATE + p % time = p % timeMax + exit DTLoop + + else ! Dist == dColl + ! Check for real or virtual collision + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + if (p % pRNG % get() < sigmaT * majorant_inv) exit DTLoop + ! Update grid distance + dGrid = dGrid - dColl + + end if + + ! Switch to surface tracking if needed + if (sigmaT * majorant_inv < ONE - self % cutoff) then + call self % surfaceTracking(p) + return + end if + + + end do DTLoop + + end subroutine deltaTracking + + !! + !! Return the inverse majorant opacity + !! For DT or HT this will be constant, for GT this will be dependent on position + !! + !! Args: + !! p [in] -> particle + !! + !! Result: + !! maj_inv -> 1 / majorant opacity + !! + function getMajInv(self, p) result (maj_inv) + class(transportOperatorTimeHT), intent(in) :: self + class(particle), intent(in) :: p + real(defReal) :: maj_inv + + if (self % method == GT) then + maj_inv = ONE / self % grid % getValue(p % coords % lvl(1) % r, p % coords % lvl(1) % dir) + else + maj_inv = ONE / self % xsData % getMajorantXS(p) + end if + + end function getMajInv + + !! + !! Provide transport operator with delta tracking/surface tracking cutoff + !! + !! Cutoff of 1 gives exclusively delta tracking, cutoff of 0 gives exclusively surface tracking + !! + subroutine init(self, dict) + class(transportOperatorTimeHT), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(nameLen) :: method + class(dictionary),pointer :: tempdict + character(100), parameter :: Here = "init (transportOperatorTimeHT_class.f90)" + + ! Initialise superclass + call init_super(self, dict) + + ! Get tracking method + call dict % getOrDefault(method, 'method', 'HT') + + select case (method) + + ! Hybrid tracking + case ('HT') + self % method = HT + ! Get cutoff value + call dict % get(self % cutoff, 'cutoff') + + ! Grid tracking + case ('GT') + self % method = GT + ! Get cutoff value + call dict % get(self % cutoff, 'cutoff') + + ! Initialise grid for hybrid tracking + tempDict => dict % getDictPtr('grid') + allocate(self % grid) + call self % grid % init(tempDict) + + ! Surface tracking + case ('ST') + self % method = ST + + ! Delta tracking + case ('DT') + self % method = DT + self % cutoff = ONE + + case default + call fatalError(Here, 'Invalid tracking method given. Must be HT, ST or DT.') + + end select + + end subroutine init + +end module transportOperatorTimeHT_class diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index 8d0e30e68..b47712343 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -19,6 +19,8 @@ module transportOperator_inter use nuclearDataReg_mod, only : ndReg_get => get use nuclearDatabase_inter, only : nuclearDatabase + ! Geometry interfaces + use trackingGrid_class, only : trackingGrid implicit none @@ -48,6 +50,8 @@ module transportOperator_inter !! Geometry pointer -> public so it can be used by subclasses (protected member) class(geometry), pointer :: geom => null() + class(trackingGrid), pointer :: grid => null() + contains ! Public interface procedure, non_overridable :: transport @@ -121,8 +125,8 @@ end subroutine transport !! Initialise transport operator from dictionary and geometry !! subroutine init(self, dict) - class(transportOperator), intent(inout) :: self - class(dictionary), intent(in) :: dict + class(transportOperator), intent(inout) :: self + class(dictionary), intent(in) :: dict ! Do nothing