diff --git a/CMakeLists.txt b/CMakeLists.txt index 61d595b91..08ed44201 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -80,6 +80,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..7014e5f4b --- /dev/null +++ b/CollisionOperator/CollisionProcessors/IMCMGstd_class.f90 @@ -0,0 +1,227 @@ +module IMCMGstd_class + + use numPrecision + use endfConstants + use universalVariables, only : VOID_MAT + 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, P_MATERIAL + 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 or MATERIAL 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") + + ! Confirm that particle is in valid material + if (p % matIdx() == VOID_MAT) call fatalError(Here, 'Collision in void material') + + ! 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) + + call p % point(dir) + + ! Sample new frequency + p % G = self % xsData % sampleEnergyGroup(p % matIdx(), p % pRNG) + + 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 % type = P_MATERIAL + + 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 83692f38f..79fa4f0c8 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 @@ -23,7 +24,8 @@ module collisionProcessorFactory_func ! For now it is necessary to adjust trailing blanks so all enteries have the same length character(nameLen),dimension(*),parameter :: AVALIBLE_collisionProcessors = [ 'neutronCEstd',& 'neutronCEimp',& - 'neutronMGstd'] + 'neutronMGstd',& + 'IMCMGstd '] contains @@ -54,6 +56,9 @@ subroutine new_collisionProcessor(new,dict) case('neutronMGstd') allocate(neutronMGstd :: new) + case('IMCMGstd') + allocate(IMCMGstd :: new) + case default print *, AVALIBLE_collisionProcessors call fatalError(Here, 'Unrecognised type of collisionProcessor: ' // trim(type)) diff --git a/CollisionOperator/collisionOperator_class.f90 b/CollisionOperator/collisionOperator_class.f90 index 22cfa1b68..4fea6bb79 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 d403b6731..6acbc393e 100644 --- a/Geometry/CMakeLists.txt +++ b/Geometry/CMakeLists.txt @@ -9,6 +9,8 @@ add_sources( ./csg_class.f90 ./geometry_inter.f90 ./geometryStd_class.f90 ./geometryReg_mod.f90 + ./discretiseGeom_class.f90 + ./geometryGrid_class.f90 ./geometryFactory_func.f90 ./fieldFactory_func.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 ca950e9cd..d9d86dd3a 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 @@ -342,6 +343,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..1de86bf70 --- /dev/null +++ b/Geometry/discretiseGeom_class.f90 @@ -0,0 +1,303 @@ + + + +!! +!! NOTE: With the addition of geometryGrid_class this module is obsolete. I left it in in case any +!! parts of it might have other uses but feel free to delete if not +!! + + + +!! +!! 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_geomIdx => geomIdx, & + gr_kill => kill + use geometryFactory_func, only : new_geometry + + ! 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 new_geometry(tempDict, geomName) + 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)//';}' + + ! Material + write(22, '(8A)') 'm'//numToChar(i)//' {temp '//numToChar(mm_matTemp(matIdx))//'; & + &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/geometryFactory_func.f90 b/Geometry/geometryFactory_func.f90 index 599f5ccf8..4b02dd07e 100644 --- a/Geometry/geometryFactory_func.f90 +++ b/Geometry/geometryFactory_func.f90 @@ -4,24 +4,26 @@ module geometryFactory_func use numPrecision - use genericProcedures, only : fatalError - use dictionary_class, only : dictionary - use charMap_class, only : charMap + use genericProcedures, only : fatalError + use dictionary_class, only : dictionary + use charMap_class, only : charMap ! Geometry - use geometry_inter, only : geometry - use geometryStd_class, only : geometryStd - use geometryReg_mod, only : gr_addGeom => addGeom + use geometry_inter, only : geometry + use geometryStd_class, only : geometryStd + use geometryGrid_class, only : geometryGrid + use geometryReg_mod, only : gr_addGeom => addGeom ! Meterial interface - use materialMenu_mod, only : mm_nameMap => nameMap + use materialMenu_mod, only : mm_nameMap => nameMap implicit none private !! Parameters - character(nameLen), dimension(*), parameter :: AVAILABLE_GEOMETRIES = ['geometryStd'] + character(nameLen), dimension(*), parameter :: AVAILABLE_GEOMETRIES = ['geometryStd ' ,& + 'geometryGrid'] ! Public interface public :: new_geometry @@ -65,6 +67,9 @@ subroutine new_geometry(dict, name, silent) case ('geometryStd') allocate(geometryStd :: geom) + case ('geometryGrid') + allocate(geometryGrid :: geom) + case default print '(A)', 'AVAILABLE GEOMETRIES' print '(A)', AVAILABLE_GEOMETRIES diff --git a/Geometry/geometryGrid_class.f90 b/Geometry/geometryGrid_class.f90 new file mode 100644 index 000000000..c5a3e0292 --- /dev/null +++ b/Geometry/geometryGrid_class.f90 @@ -0,0 +1,705 @@ +module geometryGrid_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, numToChar + use coord_class, only : coordList, coord + use dictionary_class, only : dictionary + use charMap_class, only : charMap + 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 + use materialMenu_mod, only : mm_matTemp => matTemp ,& + mm_matFile => matFile ,& + mm_init => init ,& + mm_kill => kill,& + mm_matName => matName + + implicit none + private + + !! + !! A very simplified geometry model consisting of a uniform grid. This simplicity means we don't + !! need to worry about nesting levels or complex shapes making things in general much faster. + !! Useful for many IMC benchmarks when we need to split each material region into discrete zones + !! to accurately model a time-evolving temperature field, especially when using a large number of + !! zones. + !! + !! Makes use of csg_class and copies of geometryStd_class functions for initialisation only, such + !! that input files can be virtually identical to geometryStd. + !! + !! Sample Dictionary Input: + !! geometry { + !! type geometryGrid; + !! dimensions (10 10 10); + !! + !! } + !! + !! This sample input will build the CSG geometry as in geometryStd_class, then construct a simple + !! grid geometry automatically by with 10x10x10 cells, with the material in each cell equal to the + !! material at the central point of that grid cell in the original CSG geometry. Each grid cell + !! has a new instance of each material, even if multiple grid cells are contained in a single CSG + !! material. MatIdxs will line up with localID if there are no VOID regions, if void is present + !! then this will become out of sync - however matName e.g. 'mat63' correctly corresponds to + !! position 63 (can be converted to ijk position similar to get_ijk in latUniverse_class). + !! + !! Interface: + !! Geometry Interface + !! + type, public, extends(geometry) :: geometryGrid + type(csg) :: geom + integer(shortInt), dimension(:), allocatable :: latSizeN + real(defReal), dimension(3) :: latPitch + real(defReal), dimension(3) :: corner + integer(shortInt), dimension(:,:,:), allocatable :: mats + real(defReal), dimension(6) :: geomBounds + integer(shortInt), dimension(:), allocatable :: boundary + + contains + ! Superclass procedures + procedure :: init + procedure :: kill + procedure :: placeCoord + procedure :: whatIsAt + procedure :: bounds + procedure :: move_noCache + procedure :: move_withCache + procedure :: moveGlobal + procedure :: teleport + procedure :: activeMats + + ! Public procedure unique to this class + procedure :: matBounds + + ! Private procedures + procedure, private :: initGridOnly + procedure, private :: explicitBC + procedure, private :: csg_diveToMat + procedure, private :: csg_placeCoord + procedure, private :: csg_whatIsAt + end type geometryGrid + +contains + + !! + !! Initialise geometry + !! + !! See geometry_inter for details + !! + subroutine init(self, dict, mats, silent) + class(geometryGrid), intent(inout) :: self + class(dictionary), intent(in) :: dict + type(charMap), intent(in) :: mats + logical(defBool), optional, intent(in) :: silent + logical(defBool) :: loud + real(defReal), dimension(6) :: bounds + class(surface), pointer :: surf + real(defReal), dimension(3) :: r + type(dictionary) :: matDict, tempDict + real(defReal) :: volume + character(nameLen) :: gridOnly + integer(shortInt) :: i, j, k, z, N, matIdx, uniqueID, idxCounter, voidCounter + character(100), parameter :: Here = 'init (geometryGrid_class.f90)' + + ! Choose whether to display messages + if (present(silent)) then + loud = .not.silent + else + loud = .true. + end if + + ! Get geometry discretisation + call dict % get(self % latSizeN, 'dimensions') + if (size(self % latSizeN) /= 3) call fatalError(Here, 'Dimensions must be of size 3') + + ! Allocate space for material indexes + allocate(self % mats(self % latSizeN(1),self % latSizeN(2),self % latSizeN(3))) + + ! Get boundary conditions + call dict % get(self % boundary, 'boundary') + if (size(self % boundary) /= 6) call fatalError(Here, 'boundary should be an array of size 6') + + ! Determine whether to completely rebuild geometry or to just build a grid + call dict % getOrDefault(gridOnly, 'gridOnly','n') + if (gridOnly == 'y') then + ! Switch to dedicated subroutine + call self % initGridOnly(dict) + return + else if (gridOnly /= 'n') then + call fatalError(Here, 'Unrecognised value for gridOnly setting. Should be y or n.') + end if + + ! Build the representation using CSG geometry + call self % geom % init(dict, mats, silent) + + if (loud) then + print *, "/\/\ CONVERTING CSG GEOMETRY TO GRID GEOMETRY /\/\" + end if + + ! Get geometry bounds + surf => self % geom % surfs % getPtr(self % geom % borderIdx) + bounds = surf % boundingBox() + self % geomBounds = bounds + + do i = 1, 3 + self % latPitch(i) = (bounds(i+3) - bounds(i)) / self % latSizeN(i) + end do + self % corner = bounds(1:3) + volume = product(self % latPitch) + + ! Initialise dictionary of materials for materialMenu_mod initialisation + N = product(self % latSizeN) + call matDict % init(1) + + ! Loop through each grid cell + idxCounter = 0 + voidCounter = 0 + do k = 1, self % latSizeN(3) + + ! Flip in z axis for consistency with latUniverse_class + z = self % latSizeN(3) - k + 1 + + do j = 1, self % latSizeN(2) + do i = 1, self % latSizeN(1) + + ! Get material at cell centre + r = self % corner + [i-HALF,j-HALF,z-HALF]*self % latPitch + call self % csg_whatIsAt(matIdx, uniqueID, r) + + ! Don't create a new material for void regions + if (matIdx == VOID_MAT) then + self % mats(i,j,z) = VOID_MAT + ! Count number of void regions so that matName will match up with lattice position + ! => easier data processing after obtaining results + voidCounter = voidCounter + 1 + cycle + end if + + ! Next matIdx + idxCounter = idxCounter + 1 + + ! Store in dictionary of new materials + call tempDict % init(3) + call tempDict % store('temp', mm_matTemp(matIdx)) + call tempDict % store('volume', volume) + call tempDict % store('xsFile', mm_matFile(matIdx)) + call matDict % store('mat'//numToChar(idxCounter+voidCounter), tempDict) + call tempDict % kill() + + ! Store matIdx in material array + self % mats(i,j,z) = idxCounter + + end do + end do + end do + + ! Kill current geometry and materials + call self % geom % kill() + call mm_kill() + + ! Initialise new materials + call mm_init(matDict) + + ! Kill material dictionary + call matDict % kill() + + ! Print finish line + if (loud) then + print *, "\/\/ FINISHED BUILDING GRID GEOMETRY \/\/" + print *, repeat('<>', MAX_COL/2) + end if + + end subroutine init + + !! + !! Generate grid geometry without giving any consideration to materials as is done in the full + !! subroutine above. Assigns a unique value to each grid cell so that self % whatIsAt(matIdx) + !! will return useful value. + !! + !! Unlike in full subroutine above, "bounds" is required in input dictionary. + !! + subroutine initGridOnly(self, dict) + class(geometryGrid), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal), dimension(:), allocatable :: bounds + integer(shortInt) :: i, j, k, idxCounter + + ! Get grid bounds and calculate basic properties + call dict % get(bounds, 'bounds') + self % geomBounds = bounds + do i = 1, 3 + self % latPitch(i) = (bounds(i+3) - bounds(i)) / self % latSizeN(i) + end do + self % corner = bounds(1:3) + + ! Assign a unique value to each grid cell + idxCounter = 1 + do k = 1, self % latSizeN(3) + do j = 1, self % latSizeN(2) + do i = 1, self % latSizeN(1) + self % mats(i,j,k) = idxCounter + idxCounter = idxCounter + 1 + end do + end do + end do + + end subroutine initGridOnly + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(geometryGrid), intent(inout) :: self + + call self % geom % kill() + deallocate(self % mats) + + end subroutine kill + + !! + !! Place coordinate list into geometry + !! + !! See geometry_inter for details + !! + subroutine placeCoord(self, coords) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + integer(shortInt) :: matIdx, uniqueID + character(100), parameter :: Here = 'placeCoord (geometryGrid_class.f90)' + + call self % whatIsAt(matIdx, uniqueID, coords % lvl(1) % r) + coords % matIdx = matIdx + + ! Extra unnecessary info for coords % isPlaced to return true + coords % uniqueID = matIdx + coords % nesting = 1 + + end subroutine placeCoord + + !! + !! Find material and unique cell at a given location + !! Optional direction input is not used to nudge particles across cells. All moves in this class + !! are increased by NUDGE so particles should never be exactly on a surface. + !! + !! See geometry_inter for details + !! + subroutine whatIsAt(self, matIdx, uniqueID, r, u) + class(geometryGrid), intent(in) :: self + integer(shortInt), intent(out) :: matIdx + integer(shortInt), intent(out) :: uniqueID + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), optional, intent(in) :: u + integer(shortInt), dimension(3) :: ijk + + ! Avoid rounding error for leaked particles very close to boundary + if (any(r < self % corner) .or. any(r > self % geomBounds(4:6))) then + matIdx = OUTSIDE_MAT + uniqueID = matIdx + return + end if + + ! Determine ijk location of cell + ijk = floor((r - self % corner) / self % latPitch) + 1 + + ! Get matIdx from array + matIdx = self % mats(ijk(1),ijk(2),ijk(3)) + + ! UniqueID not needed, set equal to matIdx to avoid compiler warning + uniqueID = matIdx + + end subroutine whatIsAt + + !! + !! Return Axis Aligned Bounding Box encompassing the geometry + !! + !! See geometry_inter for details + !! + function bounds(self) + class(geometryGrid), intent(in) :: self + real(defReal), dimension(6) :: bounds + + bounds = self % geomBounds + + end function bounds + + !! + !! Given coordinates placed in the geometry move point through the geometry + !! + !! See geometry_inter for details + !! + !! Uses explicit BC + !! + subroutine move_noCache(self, coords, maxDist, event) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + real(defReal) :: dist + real(defReal), dimension(3) :: r, u, r_bar + character(100), parameter :: Here = 'move (geometryGrid_class.f90)' + + ! Calculate distance to next cell crossing + r = coords % lvl(1) % r + u = coords % lvl(1) % dir + r_bar = (r - self % corner) / self % latPitch ! Normalise position within grid + r_bar = r_bar - floor(r_bar) ! Normalise position within cell + r_bar = (HALF - r_bar + sign(HALF, u)) * self % latPitch + dist = minval(r_bar / u) ! Which direction will result in crossing + + if (maxDist < dist) then ! Moves within cell + call coords % moveGlobal(maxDist) + event = COLL_EV + maxDist = maxDist ! Left for explicitness and compiler + + ! Place coords back into geometry + call self % placeCoord(coords) + + else ! Move to next cell, increased by NUDGE to avoid numerical issues + call coords % moveGlobal(dist + NUDGE) + maxDist = dist + NUDGE + + ! Set matIdx + call self % placeCoord(coords) + + ! Apply boundary conditions if leaving geometry + if (coords % matIdx == OUTSIDE_MAT) then + event = BOUNDARY_EV + call self % explicitBC(coords) + + else + ! Cell crossing within geometry - no BCs needed + event = CROSS_EV + + end if + + end if + + end subroutine move_noCache + + !! + !! Given coordinates placed in the geometry move point through the geometry + !! + !! See geometry_inter for details + !! + !! Uses explicit BC + !! + subroutine move_withCache(self, coords, maxDist, event, cache) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + type(distCache), intent(inout) :: cache + character(100), parameter :: Here = 'move_withCache (geometryGrid_class.f90)' + + ! Unnecessary + call fatalError(Here, 'Should not be called') + + end subroutine move_withCache + + !! + !! Move a particle in the top (global) level in the geometry + !! + !! See geometry_inter for details + !! + !! Uses explicit BC + !! + subroutine moveGlobal(self, coords, maxDist, event) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + real(defReal) :: dist + real(defReal), dimension(3) :: r, u, r_bar, geomSize + + ! Calculate distance to next cell crossing + r = coords % lvl(1) % r + u = coords % lvl(1) % dir + geomSize = self % geomBounds(4:6) - self % geomBounds(1:3) + r_bar = -r + self % corner + (HALF + sign(HALF, u)) * geomSize + dist = minval(r_bar / u) + + if (maxDist < dist) then ! Moves within geometry + call coords % moveGlobal(maxDist) + event = COLL_EV + maxDist = maxDist ! Left for explicitness and compiler + + ! Place coords back into geometry + call self % placeCoord(coords) + + else ! Hit geometry bounds, increased by NUDGE to avoid numerical issues + call coords % moveGlobal(dist + NUDGE) + event = BOUNDARY_EV + maxDist = dist + NUDGE + + ! Apply boundary conditions + call self % explicitBC(coords) + + end if + + + end subroutine moveGlobal + + !! + !! Move a particle in the top level without stopping + !! + !! See geometry_inter for details + !! + !! Uses co-ordinate transform boundary XSs + !! + subroutine teleport(self, coords, dist) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(in) :: dist + + ! Move the coords above the geometry + call coords % moveGlobal(dist) + + ! Place coordinates back into geometry + call self % placeCoord(coords) + + ! If point is outside apply boundary transformations + if (coords % matIdx == OUTSIDE_MAT) then + call self % explicitBC(coords) + end if + + end subroutine teleport + + !! + !! CSG equivalent normally in surface class, since we are not using surfaces perform boundary + !! transformations here instead + !! + !! Note indexing difference between bounds and boundary: + !! bounds = [xmin, ymin, zmin, xmax, ...] + !! boundary = [x-, x+, y-, y+, z-,z+] + !! + subroutine explicitBC(self, coords) + class(geometryGrid), intent(in) :: self + class(coordList), intent(inout) :: coords + integer(shortInt) :: i + real(defReal) :: outside, move + + ! Loop through axes + do i = 1, 3 + move = ZERO + + ! Negative side of bounding box + outside = (self % geomBounds(i) - coords % lvl(1) % r(i)) + if (outside >= ZERO .and. self % boundary(2*i-1) == 1) move = outside + + ! Positive side of bounding box + outside = (coords % lvl(1) % r(i) - self % geomBounds(i+3)) + if (outside >= ZERO .and. self % boundary(2*i) == 1) move = outside + + ! Move if necessary + if (move > ZERO) then + ! Flip direction + coords % lvl(1) % dir(i) = -coords % lvl(1) % dir(i) + ! Move back into geometry + call self % teleport(coords, 2*move/abs(coords % lvl(1) % dir(i))) + end if + + end do + + end subroutine explicitBC + + !! + !! Returns the list of active materials used in the geometry + !! + !! See geometry_inter for details + !! + !! NOTE: This function uses VOID_MAT and UNDEF_MAT from universalVariables + !! VOID_MAT will be listed multiple times if it occurs in multiple locations + !! + function activeMats(self) result(matList) + class(geometryGrid), intent(in) :: self + integer(shortInt), dimension(:), allocatable :: matList + + !TODO: For some reason this gives a warning after compiling, can't figure out why? + matList = reshape(self % mats,[size(self % mats)]) + + end function activeMats + + !! + !! Return position bounds of a material matIdx + !! + !! Result: + !! matBounds -> [xmin,ymin,zmin,xmax,ymax,zmax] + !! + function matBounds(self, matIdx) + class(geometryGrid), intent(in) :: self + integer(shortInt), intent(in) :: matIdx + real(defReal), dimension(6) :: matBounds + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: i, localID, temp, base + character(nameLen) :: matName + character(100), parameter :: Here = 'matBounds (geometryGrid_class.f90)' + + ! Convert matIdx to positional ID using name (will be different if void regions exist) + matName = mm_matName(matIdx) + read (matName(4:), '(I10)') localID + + ! Convert positional ID to ijk position (same as in latUniverse_class) + temp = localID - 1 + base = temp / self % latSizeN(1) + ijk(1) = temp - self % latSizeN(1) * base + 1 + temp = base + base = temp / self % latSizeN(2) + ijk(2) = temp - self % latSizeN(2) * base + 1 + ijk(3) = base + 1 + + ! Confirm that position is correct + if (self % mats(ijk(1),ijk(2),ijk(3)) /= matIdx) then + call fatalError(Here, 'Obtained matIdx different to requested matIdx') + end if + + ! Set bounds of material + do i=1, 3 + matBounds(i) = (ijk(i)-1) * self % latPitch(i) + self % geomBounds(i) + SURF_TOL + matBounds(i+3) = ijk(i) * self % latPitch(i) + self % geomBounds(i) - SURF_TOL + end do + + end function matBounds + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! CSG procedures used for initialisation +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Descend down the geometry structure untill material is reached + !! + !! Requires strting level to be specified. + !! It is private procedure common to all movment types in geometry. + !! + !! Args: + !! coords [inout] -> CoordList of a particle. Assume thet coords are already valid for all + !! levels above and including start + !! start [in] -> Starting level for meterial search + !! + !! Errors: + !! fatalError if material cell is not found untill maximum nesting is reached + !! + subroutine csg_diveToMat(self, coords, start) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + integer(shortInt), intent(in) :: start + integer(shortInt) :: rootID, localID, fill, id, i + class(universe), pointer :: uni + real(defReal), dimension(3) :: offset + character(100), parameter :: Here = 'diveToMat (geometryGrid_class.f90)' + + do i = start, HARDCODED_MAX_NEST + ! Find cell fill + rootId = coords % lvl(i) % uniRootID + localID = coords % lvl(i) % localID + call self % geom % graph % getFill(fill, id, rootID, localID) + + if (fill >= 0) then ! Found material cell + coords % matIdx = fill + coords % uniqueID = id + return + + else ! Universe fill descend a level + if (i == HARDCODED_MAX_NEST) exit ! If there is nested universe at the lowest level + + fill = abs(fill) + + ! Get current universe + uni => self % geom % unis % getPtr_fast(coords % lvl(i) % uniIdx) + + ! Get cell offset + offset = uni % cellOffset(coords % lvl(i)) + + ! Get nested universe + uni => self % geom % unis % getPtr_fast(fill) + + ! Enter nested univers + call coords % addLevel() + call uni % enter(coords % lvl(i+1), coords % lvl(i) % r - offset, coords % lvl(i) % dir) + coords % lvl(i+1) % uniRootID = id ! Must be after enter where coord has intent out + + end if + end do + + call fatalError(Here, 'Failed to find material cell. Should not happen after & + &geometry checks during build...') + + end subroutine csg_diveToMat + + + !! + !! Place coordinate list into geometry + !! + !! See geometry_inter for details + !! + subroutine csg_placeCoord(self, coords) + class(geometryGrid), intent(in) :: self + type(coordList), intent(inout) :: coords + class(universe), pointer :: uni + real(defReal), dimension(3) :: r, dir + character(100), parameter :: Here = 'placeCoord (geometryGrid_class.f90)' + + ! Check that coordList is initialised + if (coords % nesting < 1) then + call fatalError(Here, 'CoordList is not initialised. Nesting is: '//& + numToChar(coords % nesting)) + end if + + ! Place coordinates above geometry (in case they were placed) + call coords % takeAboveGeom() + + ! Enter root universe + r = coords % lvl(1) % r + dir = coords % lvl(1) % dir + + uni => self % geom % unis % getPtr_fast(self % geom % rootIdx) + + call uni % enter(coords % lvl(1), r, dir) + + coords % lvl(1) % uniRootID = 1 + + ! Dive to material + call self % csg_diveToMat(coords, 1) + + end subroutine csg_placeCoord + + + !! + !! Find material and unique cell at a given location + !! + !! See geometry_inter for details + !! + subroutine csg_whatIsAt(self, matIdx, uniqueID, r, u) + class(geometryGrid), intent(in) :: self + integer(shortInt), intent(out) :: matIdx + integer(shortInt), intent(out) :: uniqueID + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), optional, intent(in) :: u + type(coordList) :: coords + real(defReal), dimension(3) :: u_l + + ! Select direction + if (present(u)) then + u_l = u + else + u_l = [ONE, ZERO, ZERO] + end if + + ! Initialise coordinates + call coords % init(r, u_l) + + ! Place coordinates + call self % csg_placeCoord(coords) + + ! Return material & uniqueID + matIdx = coords % matIdx + uniqueID = coords % uniqueID + + end subroutine csg_whatIsAt + + +end module geometryGrid_class diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index c5ba6e00d..f7bd1cefe 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 diff --git a/InputFiles/IMC/DataFiles/hohlraumData b/InputFiles/IMC/DataFiles/hohlraumData new file mode 100644 index 000000000..8b349c73f --- /dev/null +++ b/InputFiles/IMC/DataFiles/hohlraumData @@ -0,0 +1,3 @@ + +equations hohlraum; + diff --git a/InputFiles/IMC/DataFiles/marshakData b/InputFiles/IMC/DataFiles/marshakData new file mode 100644 index 000000000..766d2a8d3 --- /dev/null +++ b/InputFiles/IMC/DataFiles/marshakData @@ -0,0 +1,3 @@ + +equations marshak; + diff --git a/InputFiles/IMC/DataFiles/olsonData b/InputFiles/IMC/DataFiles/olsonData new file mode 100644 index 000000000..e3dc09555 --- /dev/null +++ b/InputFiles/IMC/DataFiles/olsonData @@ -0,0 +1,7 @@ +// Data for 1D Olson benchamrk problem +// + +equations olson1D; + +tol 8; + diff --git a/InputFiles/IMC/DataFiles/sampleData b/InputFiles/IMC/DataFiles/sampleData new file mode 100644 index 000000000..e212109a7 --- /dev/null +++ b/InputFiles/IMC/DataFiles/sampleData @@ -0,0 +1,23 @@ + +// +// Sample material data file for IMC calculations +// Very simple, only exists to keep material creation more consistent with other calculation types +// + + + // Set of equations to use for opacities and heat capacities, currently have densmore, marshak and + // holraum. Adding new cases requires recompiling but this allows for complex multi-frequency eqns. +equations densmore; + + +/\/\ Optional Inputs /\/\ + + // Multiple either opacity or heat capacity by a constant. Useful if wanting to change constant + // value in complex equations without having to recompile (e.g. different densmore MF cases) + // Default to 1 +sigmaMultiple 100; +cvMultiple 1; + + // Constant in Fleck factor, defaults to 1 +alpha = 1; + diff --git a/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataMid b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataMid new file mode 100644 index 000000000..64badbc1d --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataMid @@ -0,0 +1,7 @@ +// Data for 1D Olson benchamrk problem +// + +equations densmore; + +sigmaMultiple 100; + diff --git a/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThick b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThick new file mode 100644 index 000000000..5fecc1d40 --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThick @@ -0,0 +1,7 @@ +// Data for 1D Olson benchamrk problem +// + +equations densmore; + +sigmaMultiple 1000; + diff --git a/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThin b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThin new file mode 100644 index 000000000..37e577cb2 --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/DataFiles/densmoreDataThin @@ -0,0 +1,7 @@ +// Data for 1D Olson benchamrk problem +// + +equations densmore; + +sigmaMultiple 10; + diff --git a/InputFiles/IMC/DensmoreMF/densmoreMid b/InputFiles/IMC/DensmoreMF/densmoreMid new file mode 100644 index 000000000..66363399c --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/densmoreMid @@ -0,0 +1,82 @@ + +type implicitPhysicsPackage; + +method IMC; +pop 100000; +limit 300000; +steps 100; +timeStep 0.01; +units ns; +printUpdates 5; + +energyGrid { + grid log; + size 30; + min 0.0001; + max 100; + } + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTime; + + // + // When running with transportOperatorGeomHT with the settings below, occasionally get error + // 'segmentation fault (core dumped)' with no further information, unable to figure out why. + // + //type transportOperatorGeomHT; + cutoff 0.9; + geometry { type geometryGrid; boundary (0 1 1 1 1 1); dimensions (25 1 1); + bounds (-2.5 -0.5 -0.5 2.5 0.5 0.5); gridOnly y; } + searchN (250 1 1); + } + +matSource { type materialSource; method fast; } + + +source { type blackBodySource; distribution surface; surface -x; temp 1; } + + +tally { + } + +geometry { + type geometryGrid; + dimensions (250 1 1); + boundary (0 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2.5 0.5 0.5); } + } + + cells {} + + universes + { + root { id 100; type rootUniverse; border 1; fill mat1; } + } + +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.001; composition {} xsFile ./DataFiles/densmoreDataMid; } + + } + +} + + + diff --git a/InputFiles/IMC/DensmoreMF/densmoreThick b/InputFiles/IMC/DensmoreMF/densmoreThick new file mode 100644 index 000000000..8e98aa253 --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/densmoreThick @@ -0,0 +1,72 @@ + +type implicitPhysicsPackage; + +method IMC; +pop 100000; +limit 300000; +steps 100; +timeStep 0.01; +units ns; +printUpdates 0; + +energyGrid { + grid log; + size 30; + min 0.0001; + max 100; + } + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTime; + } + +matSource { type materialSource; method fast; } + + +source { type blackBodySource; distribution surface; surface -x; temp 1; } + + +tally { + } + +geometry { + type geometryGrid; + dimensions (250 1 1); + boundary (0 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2.5 0.5 0.5); } + } + + cells {} + + universes + { + root { id 100; type rootUniverse; border 1; fill mat1; } + } + +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.001; composition {} xsFile ./DataFiles/densmoreDataThick; } + + } + +} + + + diff --git a/InputFiles/IMC/DensmoreMF/densmoreThin b/InputFiles/IMC/DensmoreMF/densmoreThin new file mode 100644 index 000000000..14c52807c --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/densmoreThin @@ -0,0 +1,72 @@ + +type implicitPhysicsPackage; + +method IMC; +pop 100000; +limit 300000; +steps 100; +timeStep 0.01; +units ns; +printUpdates 0; + +energyGrid { + grid log; + size 30; + min 0.0001; + max 100; + } + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTime; + } + +matSource { type materialSource; method fast; } + + +source { type blackBodySource; distribution surface; surface -x; temp 1; } + + +tally { + } + +geometry { + type geometryGrid; + dimensions (250 1 1); + boundary (0 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2.5 0.5 0.5); } + } + + cells {} + + universes + { + root { id 100; type rootUniverse; border 1; fill mat1; } + } + +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.001; composition {} xsFile ./DataFiles/densmoreDataThin; } + + } + +} + + + diff --git a/InputFiles/IMC/DensmoreMF/twoRegion1 b/InputFiles/IMC/DensmoreMF/twoRegion1 new file mode 100644 index 000000000..91902af59 --- /dev/null +++ b/InputFiles/IMC/DensmoreMF/twoRegion1 @@ -0,0 +1,78 @@ + +type implicitPhysicsPackage; + +method IMC; +pop 10000; +limit 30000; +steps 100; +timeStep 0.01; +units ns; +printUpdates 0; + +energyGrid { + grid log; + size 30; + min 0.0001; + max 100; + } + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTime; + } + +matSource { type materialSource; method fast; } + + +source { type blackBodySource; distribution surface; surface -x; temp 1; } + + +tally { + } + +geometry { + type geometryGrid; + dimensions (250 1 1); + boundary (0 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 1.5 0.5 0.5); } + sep { id 2; type xPlane; x0 0.5; } + } + + cells { + call1 {id 1; type simpleCell; surfaces (-2); filltype mat; material mat1; } + call2 {id 2; type simpleCell; surfaces (2); filltype mat; material mat2; } + } + + universes + { + root { id 1; type rootUniverse; border 1; fill u<2>; } + cells { id 2; type cellUniverse; cells (1 2); } + } + +} + +nuclearData { + + handles { + mg { type baseMgIMCDatabase; } + } + + + materials { + + mat1 { temp 0.001; composition {} xsFile ./DataFiles/densmoreDataThin; } + mat2 { temp 0.001; composition {} xsFile ./DataFiles/densmoreDataThick; } + + } + +} + + + diff --git a/InputFiles/IMC/hohlraum b/InputFiles/IMC/hohlraum new file mode 100644 index 000000000..8b1701bc0 --- /dev/null +++ b/InputFiles/IMC/hohlraum @@ -0,0 +1,100 @@ + +type implicitPhysicsPackage; + +method IMC; +pop 500000; +limit 2000000; +steps 100; +timeStep 0.01; +units ns; + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + //type transportOperatorTime; + type transportOperatorGeomHT; + cutoff 0.5; + geometry { type geometryGrid; boundary (0 0 0 0 1 1); dimensions (40 40 1); + bounds (-0.5 -0.5 -0.5 0.5 0.5 0.5); gridOnly y; } + searchN (200 200 1); + } + +source { + type blackBodySource; + distribution surface; + surface -x; + temp 1; +} + +viz { vizDict { type vtk; corner (-0.5 -0.5 -0.5); width (1 1 1); vox (20 20 1); } } + +tally { + } + +geometry { + type geometryGrid; + dimensions (200 200 1); + 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.000001; composition {} xsFile ./DataFiles/hohlraumData; } + + } + +} + + + diff --git a/InputFiles/IMC/marshakWave b/InputFiles/IMC/marshakWave new file mode 100644 index 000000000..216fe6912 --- /dev/null +++ b/InputFiles/IMC/marshakWave @@ -0,0 +1,42 @@ +// Marshak wave IMC benchmark +// Requires lightSpeed and radiation constant set to ONE in universalVariables + +type implicitPhysicsPackage; + +method IMC; +pop 100000; +limit 250000; +steps 1000; +timeStep 0.5; +units marshak; + +collisionOperator { photonMG {type IMCMGstd;} } + +transportOperator { + type transportOperatorTime; + } + +tally {} + +// Black body surface source +source { type blackBodySource; distribution surface; surface -x; temp 1; } + +geometry { + type geometryGrid; + dimensions (100 1 1); + 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/olson1D b/InputFiles/IMC/olson1D new file mode 100644 index 000000000..e8bd2b347 --- /dev/null +++ b/InputFiles/IMC/olson1D @@ -0,0 +1,70 @@ + +type implicitPhysicsPackage; + +method ISMC; +pop 5000; +limit 1000000; +steps 1333; +timeStep 0.0001; +units ns; +printUpdates 5; + + +energyGrid { + grid unstruct; + bins (17.8 10.0 5.62 3.16 1.78 1.0 0.562 0.316 + 0.178 0.1 0.0562 0.0316 0.0178 0.01 + 0.00562 0.00316 0.00178 0.001 0.0001); + } + +collisionOperator { + photonMG {type IMCMGstd;} + } + +transportOperator { + type transportOperatorTime; + } + +source { type blackBodySource; distribution olson1D; temp 0.5; untilStep 667; } + + +tally { + } + +geometry { + type geometryGrid; + dimensions (100 1 1); + boundary (1 1 1 1 1 1); + graph {type shrunk;} + + surfaces + { + outer { id 1; type box; origin ( 0.0 0.0 0.0); halfwidth ( 2.4 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/olsonData; } + + } + +} + + + diff --git a/InputFiles/IMC/sampleInput b/InputFiles/IMC/sampleInput new file mode 100644 index 000000000..2ffb73035 --- /dev/null +++ b/InputFiles/IMC/sampleInput @@ -0,0 +1,54 @@ + +// Intended as a sample/tutorial input file for calculations using IMC physics +// package, to detail settings that differ to other calculation types + +type implicitPhysicsPackage; + +method IMC; // IMC or ISMC +pop 10000; // +limit 25000; // +steps 5000; // Number of timesteps to take +timeStep 0.1; // Timestep for calculation +units marshak; // Units for conversion of timeStep in physicsPackage +printUpdates 2; // Prints first N brief material updates per timestep + +// Energy grid for multi-frequency calculations, if not provided then defaults to grey case +energygrid { + grid log; + size 30; + min 0.0001; + max 100; + } + +collisionOperator { photonMG {type IMCMGstd;} } + +transportOperator { + type transportOperatorTime; // Standard transport operator for IMC/ISMC + } + +tally {} // Required tallies are defined automatically in physicsPackage + +// Black body surface source placed on -x face of geometry bounding box with T = 1 +source { type blackBodySource; distribution surface; surface -x; temp 1; } + +geometry { + type geometryGrid; // Splits geometry into a uniform grid with the dimensions given below + // Efficient for IMC, see geometryGrid_class.f90 for more detail + dimensions (100 1 1); + + boundary (0 0 1 1 1 1); // 0 -> Open boundary, 1 -> Reflective boundary + 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/sampleData; }} +} + + + 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..e19d6d06e --- /dev/null +++ b/NuclearData/IMCMaterial_inter.f90 @@ -0,0 +1,105 @@ +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(getFleck), deferred :: getFleck + procedure(getTemp), deferred :: getTemp + 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 + + !! + !! 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 + + 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..45f51bd78 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -10,13 +10,15 @@ !! 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 +!! matFile -> Return file path to material data given matIdx +!! matIdx -> Return material Index given Name !! module materialMenu_mod @@ -64,6 +66,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 +96,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 +116,8 @@ module materialMenu_mod public :: getMatPtr public :: nMat public :: matName + public :: matTemp + public :: matFile public :: matIdx contains @@ -131,7 +137,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 +210,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 +227,48 @@ 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 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 !! @@ -267,13 +315,16 @@ subroutine init_materialItem(self, name, dict) ! Return to initial state call self % kill() - ! Load easy components c + ! Load easy components 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') - call compDict % keys(keys) + if (dict % isPresent('composition')) then + compDict => dict % getDictPtr('composition') + call compDict % keys(keys) + end if ! Allocate space for nuclide information allocate(self % nuclides(size(keys))) @@ -288,24 +339,24 @@ subroutine init_materialItem(self, name, dict) end if ! Load definitions - do i =1,size(keys) - ! Check if S(a,b) is on and required for that nuclide - if (hasSab .and. moderDict % isPresent(keys(i))) then - self % nuclides(i) % hasSab = .true. - call moderDict % get(self % nuclides(i) % file_Sab, keys(i)) - end if - - ! Initialise the nuclides - call compDict % get(self % dens(i), keys(i)) - call self % nuclides(i) % init(keys(i)) - end do + if (associated(compDict)) then + do i =1,size(keys) + ! Check if S(a,b) is on and required for that nuclide + if (hasSab .and. moderDict % isPresent(keys(i))) then + self % nuclides(i) % hasSab = .true. + call moderDict % get(self % nuclides(i) % file_Sab, keys(i)) + end if + + ! Initialise the nuclides + call compDict % get(self % dens(i), keys(i)) + call self % nuclides(i) % init(keys(i)) + end do + + end if ! Save dictionary self % extraInfo = dict - ! TODO: Remove composition subdictionary from extraInfo - ! Or rather do not copy it in the first place - end subroutine init_materialItem !! @@ -318,6 +369,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..2733f0557 --- /dev/null +++ b/NuclearData/mgIMCData/CMakeLists.txt @@ -0,0 +1,9 @@ +# Add source files for compilation +add_sources(./mgIMCMaterial_inter.f90 + ./mgIMCDatabase_inter.f90 + ./baseMgIMC/baseMgIMCMaterial_class.f90 + ./baseMgIMC/baseMgIMCDatabase_class.f90 + ./baseMgIMC/materialEquations.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..b7ff9391a --- /dev/null +++ b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCDatabase_class.f90 @@ -0,0 +1,526 @@ +module baseMgIMCDatabase_class + + use numPrecision + use endfConstants + use universalVariables, only : VOID_MAT, OUTSIDE_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 + use RNG_class, only : RNG + + ! 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, mm_matName => matName + + ! baseMgIMC Objects + use baseMgIMCMaterial_class, only : baseMgIMCMaterial + use materialEquations, only : mgEnergyGrid + + ! Energy grid for multi-frequency + use energyGrid_class, only : energyGrid + use energyGridRegistry_mod, only : get_energyGrid + + implicit none + private + + !! + !! Public Pointer Cast + !! + public :: baseMgIMCDatabase_TptrCast + public :: baseMgIMCDatabase_CptrCast + + !! + !! Database for MG IMC and ISMC calculations + !! + !! All materials in a problem are of class baseMgIMCMaterial. 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 + !! nG -> number of energy groups + !! + !! Interface: + !! nuclearDatabase interface + !! + type, public, extends(mgIMCDatabase) :: baseMgIMCDatabase + type(baseMgIMCMaterial), dimension(:), pointer :: mats => null() + integer(shortInt), dimension(:), allocatable :: activeMats + integer(shortInt) :: nG = 0 + type(energyGrid), private :: eGrid + + contains + ! Superclass Interface + procedure :: getTransMatXS + procedure :: getTotalMatXS + procedure :: getMajorantXS + procedure :: matNamesMap + procedure :: getMaterial + procedure :: getNuclide + procedure :: getReaction + procedure :: getEmittedRad + procedure :: getMaterialEnergy + procedure :: updateProperties + procedure :: setCalcType + procedure :: sampleTransformTime + procedure :: sampleEnergyGroup + 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 + character(100), parameter :: Here = 'getTotalMatXS (baseMgIMCDatabase_class.f90)' + + if (matIdx == OUTSIDE_MAT) call fatalError(Here, 'Requesting XS in OUTSIDE_MAT') + + 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 the 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 + + !! + !! Return material energy + !! + !! Args: + !! matIdx [in] [optional] -> If provided, return the energy of only matIdx + !! Otherwise, return the total energy of all mats + !! + function getMaterialEnergy(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) % getMatEnergy() + return + end if + + ! Otherwise, return total energy emitted from all materials + energy = 0 + + do i=1, size(self % mats) + energy = energy + self % mats(i) % getMatEnergy() + end do + + end function getMaterialEnergy + + !! + !! 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), not in parallel to allow correct order of console output + do i = 1, printUpdates + print * + print *, ' Material update for '//trim(mm_matName(i))//':' + call self % mats(i) % updateMat(tallyEnergy(i), .true.) + if (i == printUpdates) print * + 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 + + !! + !! Tell each material if we are using IMC or ISMC + !! + subroutine setCalcType(self, type) + class(baseMgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: type + integer(shortInt) :: i + + do i=1, size(self % mats) + call self % mats(i) % setCalcType(type) + end do + + end subroutine setCalcType + + !! + !! Sample energy group of a particle emitted from material matIdx + !! + !! Args: + !! matIdx [in] -> index of material to sample from + !! rand [in] -> RNG + !! + !! Result: + !! G -> energy group of sampled particle + !! + function sampleEnergyGroup(self, matIdx, rand) result(G) + class(baseMgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + class(RNG), intent(inout) :: rand + integer(shortInt) :: G + + G = self % mats(matIdx) % sampleEnergyGroup(rand) + + end function sampleEnergyGroup + + !! + !! Sample the time taken for a material particle to transform into a photon + !! Used for ISMC only + !! + function sampleTransformTime(self, matIdx, rand) result(t) + class(baseMgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + class(RNG), intent(inout) :: rand + real(defReal) :: t + + t = self % mats(matIdx) % sampleTransformTime(rand) + + end function sampleTransformTime + + !! + !! 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 + logical(defBool) :: err + character(nameLen) :: gridName + 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 + + ! Get energy grid in multi-frequency case + gridName = 'mg' + call get_energyGrid(self % eGrid, gridName, err) + if (err .eqv. .false.) then + self % nG = self % eGrid % getSize() + mgEnergyGrid => self % eGrid + else + self % nG = 1 + 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) + + ! Verify number of groups + if(self % nG /= self % mats(i) % nGroups()) then + call fatalError(Here,'Inconsistent # 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..bf2ee5521 --- /dev/null +++ b/NuclearData/mgIMCData/baseMgIMC/baseMgIMCMaterial_class.f90 @@ -0,0 +1,660 @@ +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 simulationTime_class, only : timeStep + + ! Nuclear Data Interfaces + use materialHandle_inter, only : materialHandle + use mgIMCMaterial_inter, only : mgIMCMaterial, kill_super => kill + use IMCXSPackages_class, only : IMCMacroXSs + use materialEquations + + + 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 :: EMISSION_PROB = 4 + + ! Calculation Type + integer(shortInt), parameter, public :: IMC = 1 + integer(shortInt), parameter, public :: ISMC = 2 + + !! + !! Material class for MG IMC and ISMC calculations + !! + !! Stores MG data in a table. + !! + !! Public members: + !! data -> Rank 2 array with all XSs data + !! + !! Interface: + !! materialHandle interface + !! mgIMCMaterial interface + !! init -> initialise Basic MG Material from dictionary and config keywords + !! 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 + !! getMatEnergy -> returns energy of material + !! setCalcType -> set to IMC or ISMC + !! sampleEnergyGroup -> return sampled energy group of an emitted photon + !! sampleTransformTime -> return sampled time for transform of P_MATERIAL to P_PHOTON (ISMC) + !! + !! Note: + !! Order of "data" array is: data(XS_type, Group #) + !! For multigroup calculations: + !! -> energyGrid { } is given in input file, not in material data file + !! -> if not provided, simulation is done for grey case + !! -> numberOfGroups is automatically extracted from energyGrid + !! + !! Data file sample inputs: + !! equations olson; -> required input, tells code which set of compiled equations to use + !! (see materialEquations.f90 for examples and to add new cases) + !! sigmaFactor 10; -> multiplies sigma equation by a constant, avoids recompiling or + !! duplicating equations + !! cv 1; -> as above, for heat capacity + !! tol 8; -> defaults to 6 if not given. The order of tolerence used in temp + !! calculation, i.e. 8 => 1e-8, this is then multiplied by current temp + !! alpha 0.7; -> optional parameter used to tweak fleck factor, defaults to 1 + !! + type, public, extends(mgIMCMaterial) :: baseMgIMCMaterial + real(defReal),dimension(:,:), allocatable :: data ! XS (opacity) data for each group + character(nameLen) :: name ! Name for update equations + real(defReal) :: T ! Temperature + real(defReal) :: V ! Volume + real(defReal) :: fleck ! Fleck factor + real(defReal) :: alpha ! User-defined parameter for fleck factor + real(defReal) :: sigmaP ! Planck opacity + real(defReal) :: matEnergy ! Total energy stored in material + real(defReal) :: prevMatEnergy ! Energy prior to material update + real(defReal) :: energyDens ! Energy density = matEnergy/V + real(defReal) :: eta ! aT^4/energyDens (ISMC only) + real(defReal) :: sigmaFactor ! Constant to multiply sigma by + real(defReal) :: cvFactor ! Constant to multiply heat capacity by + real(defReal) :: tol ! Tolerance for calculating temperature + integer(shortInt) :: calcType ! IMC or ISMC + + contains + ! Superclass procedures + procedure :: kill + procedure :: getMacroXSs_byG + procedure :: getTotalXS + procedure :: getFleck + procedure :: getTemp + + ! Local procedures + procedure :: init + procedure :: nGroups + procedure :: updateMat + procedure :: getEmittedRad + procedure :: getMatEnergy + procedure :: setCalcType + procedure :: sampleEnergyGroup + procedure :: sampleTransformTime + + ! Private local procedures + procedure, private :: tempFromEnergy + procedure, private :: sigmaFromTemp + procedure, private :: updateFleck + + end type baseMgIMCMaterial + +contains + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Standard procedures +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! 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) + + 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, i + real(defReal) :: dT, tempT, tempU, tol + character(100), parameter :: Here = 'init (baseMgIMCMaterial_class.f90)' + + ! Read number of groups + if (associated(mgEnergyGrid)) then + nG = mgEnergyGrid % getSize() + if (nG < 1) call fatalError(Here,'Number of groups is invalid' // numToChar(nG)) + else + nG = 1 + end if + + ! Allocate space for data + N = 4 + allocate(self % data(N, nG)) + + ! Store alpha setting + call dict % getOrDefault(self % alpha, 'alpha', ONE) + + ! Get name of equations for heat capacity and opacity calculations + call dict % get(self % name, 'equations') + + ! Get optional multiplication factor for heat capacity and opacity + call dict % getOrDefault(self % sigmaFactor, 'sigmaMultiple', ONE) + call dict % getOrDefault(self % cvFactor, 'cvMultiple', ONE) + + ! Get optional tolerance for temperature calculation + call dict % getOrDefault(tol, 'tol', 6*ONE) + self % tol = 10_defReal**(-tol) + + ! Read initial temperature and volume + call dict % get(self % T, 'T') + call dict % get(self % V, 'V') + + ! Calculate initial opacities + call self % sigmaFromTemp() + + ! Calculate initial energy density by integration + dT = self % T / 1000 + tempT = dT/2 + tempU = 0 + do i=1, 1000 + tempU = tempU + dT * evaluateCv(self % name, tempT) + tempT = tempT + dT + end do + self % energyDens = tempU * self % cvFactor + self % matEnergy = self % energyDens * self % V + + end subroutine init + + !! + !! 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 + + !! + !! Set the calculation type to be used + !! + !! Current options: + !! IMC + !! ISMC + !! + !! Errors: + !! Unrecognised option + !! + subroutine setCalcType(self, calcType) + class(baseMgIMCMaterial), intent(inout) :: self + integer(shortInt), intent(in) :: calcType + character(100), parameter :: Here = 'setCalcType (baseMgIMCMaterial_class.f90)' + + if(calcType /= IMC .and. calcType /= ISMC) call fatalError(Here, 'Invalid calculation type') + + self % calcType = calcType + + ! Set initial fleck factor + call self % updateFleck() + + end subroutine setCalcType + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Material Updates +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! 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 + !! loud [in, optional] -> Bool, if true then will print updates to screen + !! + subroutine updateMat(self, tallyEnergy, loud) + class(baseMgIMCMaterial),intent(inout) :: self + real(defReal), intent(in) :: tallyEnergy + logical(defBool), intent(in), optional :: loud + real(defReal) :: prevTemp, change + character(100), parameter :: Here = 'updateMat (baseMgIMCMaterial_class.f90)' + + ! Save previous energy and temperature + self % prevMatEnergy = self % matEnergy + prevTemp = self % T + + ! Update material internal energy + if (self % calcType == IMC) then + self % matEnergy = self % matEnergy - self % getEmittedRad() + tallyEnergy + else + self % matEnergy = tallyEnergy + end if + + self % energyDens = self % matEnergy / self % V + + ! Confirm new energy density is valid + if (self % energyDens <= ZERO) call fatalError(Here, 'Energy density is not positive') + + ! Update material temperature + self % T = self % tempFromEnergy() + + ! Update sigma + call self % sigmaFromTemp() + + ! Update fleck factor + call self % updateFleck() + + ! Print updates if requested + if (present(loud)) then + if (loud) then + change = self % matEnergy - self % prevMatEnergy + if (change < ZERO) then + print *, ' Mat Energy ='//numToChar(self % matEnergy)//& + ' ( -'//numToChar(abs(change))//' )' + else + print *, ' Mat Energy ='//numToChar(self % matEnergy)//& + ' ( +'//numToChar(change)//' )' + end if + change = self % T - prevTemp + if (change < ZERO) then + print *, ' Mat Temperature ='//numToChar(self % T)//& + ' ( -'//numToChar(abs(change))//' )' + else + print *, ' Mat Temperature ='//numToChar(self % T)//& + ' ( +'//numToChar(change)//' )' + end if + + end if + end if + + end subroutine updateMat + + !! + !! Calculate material temperature from internal energy by integration + !! + !! Step up (or down if energy is lower than in previous step) through temperature in steps of dT + !! and increment energy density by dT*cv(T) until target energy density is reached + !! + function tempFromEnergy(self) result(T) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: T, dT, tempT, U, tempU, tol, error, increase + integer(shortInt) :: i + character(100), parameter :: Here = 'tempFromEnergy (mgIMCMaterial_class.f90)' + + ! Confirm that current temperature is valid + if (self % T <= ZERO) call fatalError(Here, 'Zero temperature') + + ! Parameters affecting accuracy + dT = self % T / 1000 ! Initial temperature step size to take, reduced after overshoot + tol = self % T * self % tol ! Continue stepping until within tolerance + + ! Starting temperature and energy density + T = self % T + U = self % prevMatEnergy / self % V + + ! If starting energy density is higher than target, flip dT to be negative + if (U > self % energyDens) dT = -dT + + i = 0 + + integrate:do + + ! Protect against infinite loop + i = i+1 + if (i > 100000) then + print *, 'Energy density: ', self % energyDens + call fatalError(Here, '100,000 iterations without convergence, maybe NaN energy density?') + end if + ! Increase step size to avoid lack of convergence due to very small starting temperature + if (mod(i,1000)==0) dT = 10*dT + + ! Increment temperature and increment the corresponding energy density + tempT = T + dT/2 + increase = dT * evaluateCv(self % name, tempT) * self % cvFactor + ! Protect against division by 0 or other numerical errors + if (increase /= increase .or. increase > INF) increase = ZERO + tempU = U + increase + + error = self % energyDens - tempU + + if (abs(error) < tol) then ! Finished + T = T + dT + exit integrate + end if + + if (error*dT < 0) then ! If error and dT have different signs then we have overshot + ! Decrease dT and try again + dT = dT/2 + cycle integrate + end if + + ! Update temperature and energy and return to start + T = T + dT + U = tempU + + end do integrate + + end function tempFromEnergy + + !! + !! Calculate opacities from current temp + !! + subroutine sigmaFromTemp(self) + class(baseMgIMCMaterial), intent(inout) :: self + integer(shortInt) :: i, j, stepsPerGroup + real(defReal) :: sigmaP, E, EStep, increase, upper, lower + character(100), parameter :: Here = 'sigmaFromTemp (baseMgIMCMaterial_class.f90)' + + ! Evaluate opacities for grey case + if (self % nGroups() == 1) then + self % data(CAPTURE_XS,1) = evaluateSigma(self % name, self % T, ONE) * self % sigmaFactor + self % data(IESCATTER_XS,1) = ZERO + self % data(TOTAL_XS,1) = self % data(CAPTURE_XS,1) + self % data(IESCATTER_XS,1) + ! Planck opacity equal to absorption opacity for single frequency + self % sigmaP = self % data(CAPTURE_XS, 1) + return + end if + + ! Evaluate opacities for frequency-dependent case + do i = 1, self % nGroups() + ! Take geometric mean of upper and lower group boundaries + upper = evaluateSigma(self % name, self % T, mgEnergyGrid % bin(i)) + lower = evaluateSigma(self % name, self % T, mgEnergyGrid % bin(i+1)) + self % data(CAPTURE_XS,i) = sqrt(upper*lower)*self % sigmaFactor !min(1e5_defReal, sqrt(o1 * o2)) + + upper = upper * normPlanckSpectrum(mgEnergyGrid % bin(i), self % T) + lower = lower * normPlanckSpectrum(mgEnergyGrid % bin(i+1), self % T) + self % data(EMISSION_PROB,i) = sqrt(upper*lower) + self % data(EMISSION_PROB,i) = self % data(EMISSION_PROB,i)*(mgEnergyGrid % bin(i) - mgEnergyGrid % bin(i+1)) + + end do + self % data(IESCATTER_XS,:) = ZERO + self % data(TOTAL_XS,:) = self % data(CAPTURE_XS,:) + self % data(IESCATTER_XS,:) + self % data(EMISSION_PROB,:) = self % data(EMISSION_PROB,:)/sum(self % data(EMISSION_PROB,:)) + + ! Evaluate Planck opacity via a much finer numerical integration + sigmaP = ZERO + do i = 1, self % nGroups() + ! 100 integration steps per energy group seems to give very good accuracy, might be worth playing around with + stepsPerGroup = 100 + EStep = (mgEnergyGrid % bin(i) - mgEnergyGrid % bin(i+1)) / stepsPerGroup + E = mgEnergyGrid % bin(i) - 0.5*EStep + do j = 1, stepsPerGroup + increase = normPlanckSpectrum(E, self % T)*evaluateSigma(self % name, self % T, E) + sigmaP = sigmaP + EStep * increase + E = E - EStep + end do + end do + self % sigmaP = sigmaP * self % sigmaFactor + + end subroutine sigmaFromTemp + + !! + !! Update fleck factor + !! + subroutine updateFleck(self) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: beta, zeta + character(100), parameter :: Here = 'updateFleck (baseMgIMCMaterial_class.f90)' + + ! Calculate beta, ratio of radiation and material heat capacities + beta = 4 * radiationConstant * self % T**3 / (evaluateCv(self % name, self % T)*self % cvFactor) + + ! Use time step size to calculate fleck factor + select case(self % calcType) + + case(IMC) + self % fleck = 1/(1+self % sigmaP*lightSpeed*beta*timeStep()*self % alpha) + + case(ISMC) + self % eta = radiationConstant * self % T**4 / self % energyDens + zeta = beta - self % eta + self % fleck = 1 / (1 + zeta*self % sigmaP*lightSpeed*timeStep()) + + case default + call fatalError(Here, 'Unrecognised calculation type') + + end select + + end subroutine updateFleck + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Obtain material info +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! 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 * timeStep() * 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 + + !! + !! Return energy per unit volume of material + !! + function getMatEnergy(self) result(energy) + class(baseMgIMCMaterial), intent(inout) :: self + real(defReal) :: energy + + energy = self % matEnergy + + end function getMatEnergy + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Sample emission properties +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Sample energy group for a particle emitted from material (and after effective scattering) + !! + function sampleEnergyGroup(self, rand) result(G) + class(baseMgIMCMaterial), intent(inout) :: self + class(RNG), intent(inout) :: rand + integer(shortInt) :: G, i + real(defReal) :: random, cumProb + + if (self % nGroups() == 1) then + G = 1 + return + end if + + random = rand % get() + cumProb = ZERO + + ! Choose based on emission probability of each group + do i = 1, self % nGroups() + cumProb = cumProb + self % data(EMISSION_PROB,i) + if (random < cumProb) then + G = i + return + end if + end do + + end function sampleEnergyGroup + + !! + !! Sample the time taken for a material particle to transform into a photon + !! Used for ISMC only + !! + function sampleTransformTime(self, rand) result(t) + class(baseMgIMCMaterial), intent(inout) :: self + class(RNG), intent(inout) :: rand + real(defReal) :: t + + t = -log(rand % get()) / (self % sigmaP * self % fleck * self % eta * lightSpeed) + + end function sampleTransformTime + +end module baseMgIMCMaterial_class diff --git a/NuclearData/mgIMCData/baseMgIMC/materialEquations.f90 b/NuclearData/mgIMCData/baseMgIMC/materialEquations.f90 new file mode 100644 index 000000000..eb8f4fb32 --- /dev/null +++ b/NuclearData/mgIMCData/baseMgIMC/materialEquations.f90 @@ -0,0 +1,195 @@ +!! +!! Module to store temperature and/or frequency-dependent equations for materials, especially +!! those too complicated to be easily read in from an input file. Also contains an energy grid to +!! allow materials to access particle energy group bounds for use in evaluating these equations. +!! +!! Also stores energy grid for multigroup problems for easy access by materials +!! +!! For a new set of material equations: +!! -> Add name to AVAILABLE_equations +!! -> Add case to evaluateCv and evaluateSigma +!! -> Evaluate simple equations (e.g. 'marshak' or 'hohlraum') in these functions, +!! or can link to new functions (e.g. 'olson1D') +!! +!! Virtually all benchmark problems for IMC have a real scattering opactiy of 0, giving only +!! effective absorptions and effective scattering. If desired, incorporating real scattering should +!! be very easy, adding a new equation type in this module and scattering support in material and +!! collision operator. +!! +module materialEquations + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError + use energyGrid_class, only : energyGrid + + implicit none + private + + character(nameLen),dimension(*),parameter :: AVAILABLE_equations = ['marshak ',& + 'hohlraum',& + 'olson1D ',& + 'densmore'] + + public :: evaluateCv + public :: evaluateSigma + public :: normPlanckSpectrum + + ! Energy grid for multi-frequency problems for easy access by material classes + type(energyGrid), pointer, public :: mgEnergyGrid => null() + + interface evaluateCv + module procedure evaluateCv + end interface + + interface evaluateSigma + module procedure evaluateSigma + end interface + + + contains + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Non-specific interface equations - add new case to each +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Evaluate heat capacity equation at a given temperature + !! + !! Args: + !! equation -> name of the equation to use + !! T -> temperature to evaluate at + !! + function evaluateCv(equation, T) result(cv) + character(*), intent(in) :: equation + real(defReal), intent(in) :: T + real(defReal) :: cv + character(100), parameter :: Here = 'getCv (materialEquations.f90)' + + select case(equation) + + case('marshak') + cv = 7.14 + + case('hohlraum') + cv = 0.3 + + case('olson1D') + cv = cvOlson1D(T) + + case('densmore') + cv = 0.1 + + case default + cv = ZERO + print *, AVAILABLE_equations + call fatalError(Here, 'Unrecognised equation: '//trim(equation)) + + end select + + end function evaluateCv + + !! + !! Evaluate opacity equation at a given temperature and frequency (energy E) + !! + !! Args: + !! equation -> name of the equation to use + !! T -> temperature to evaluate at + !! nu -> frequency of the particle + !! + function evaluateSigma(equation, T, E) result(sigma) + character(*), intent(in) :: equation + real(defReal), intent(in) :: T + real(defReal), intent(in) :: E + real(defReal) :: sigma + character(100), parameter :: Here = 'getSigma (materialEquations.f90)' + + select case(equation) + + case('marshak') + sigma = 10*T**(-3) + + case('hohlraum') + sigma = 100*T**(-3) + + case('olson1D') + sigma = sigmaOlson1D(T, E) + + case('densmore') + sigma = ONE / (E**3 * sqrt(T)) + + case default + sigma = ZERO + print *, AVAILABLE_equations + call fatalError(Here, 'Unrecognised equation: '//trim(equation)) + + end select + + end function evaluateSigma + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Custom material equations for various input files +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! 1-Dimensional Multi-frequency test by Olson (2020) + !! + !! Olson, G.L., 2020. Stretched and filtered multigroup pn transport for improved positivity + !! and accuracy. Journal of Computational and Theoretical Transport 0, 1–18. + !! + function cvOlson1D(T) result(cv) + real(defReal), intent(in) :: T + real(defReal) :: cv, root, alpha, dAlphadT + + root = sqrt(1+4*exp(0.1/T)) + alpha = 0.5*exp(-0.1/T)*(root-1) + dAlphadT = 0.1*(alpha-1/root)/(T*T) + + cv = radiationConstant*0.1*(1+alpha+(T+0.1)*dAlphadT) + + ! Deal with numerical errors from poorly defined regions (e.g. T almost 0) + if (cv /= cv .or. cv > INF) cv = ZERO + + end function cvOlson1D + + function sigmaOlson1D(T, E) result(sigma) + real(defReal), intent(in) :: T + real(defReal), intent(in) :: E + real(defReal) :: sigma + + if (E < 0.008) then + sigma = min(1e7_defReal, 1e9_defReal*T*T) + else if (E < 0.3) then + sigma = 192/(E*E*(1+200*T**1.5)) + else + sigma = 192*sqrt(0.3)/(E**2.5*(1+200*T**1.5)) + 4e4_defReal*(0.3/E)**2.5/(1+8000*T**2) + end if + + ! Multiply by density + sigma = 0.001*sigma + + end function sigmaOlson1D + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Commonly used equations +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Evaluate frequency-normalised Planck spectrum + !! + !! Args: + !! nu -> frequency + !! T -> temperature + !! + pure function normPlanckSpectrum(E, T) result(b) + real(defReal), intent(in) :: E + real(defReal), intent(in) :: T + real(defReal) :: b + + b = 15*E**3 / ((pi*T)**4 * (exp(E/T)-1)) + + end function normPlanckSpectrum + + +end module materialEquations diff --git a/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 b/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 new file mode 100644 index 000000000..7d8947116 --- /dev/null +++ b/NuclearData/mgIMCData/mgIMCDatabase_inter.f90 @@ -0,0 +1,144 @@ +module mgIMCDatabase_inter + + use numPrecision + use RNG_class, only : RNG + + ! 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(getMaterialEnergy), deferred :: getMaterialEnergy + procedure(updateProperties), deferred :: updateProperties + procedure(setCalcType), deferred :: setCalcType + procedure(sampleEnergyGroup), deferred :: sampleEnergyGroup + procedure(sampleTransformTime), deferred :: sampleTransformTime + + 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 + + !! + !! Return material energy + !! + !! Args: + !! matIdx [in] [optional] -> If provided, return the energy of only matIdx + !! Otherwise, return total energy of all mats + !! + function getMaterialEnergy(self, matIdx) result(energy) + import :: mgIMCDatabase, shortInt, defReal + class(mgIMCDatabase), intent(in) :: self + integer(shortInt), intent(in), optional :: matIdx + real(defReal) :: energy + end function getMaterialEnergy + + + !! + !! 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 + + !! + !! Tell each material if we are using IMC or ISMC + !! + subroutine setCalcType(self, type) + import :: mgIMCDatabase, shortInt + class(mgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: type + end subroutine setCalcType + + !! + !! Sample energy group of a particle emitted from material matIdx + !! + !! Args: + !! matIdx [in] -> index of material to sample from + !! rand [in] -> RNG + !! + !! Result: + !! G -> energy group of sampled particle + !! + function sampleEnergyGroup(self, matIdx, rand) result(G) + import :: mgIMCDatabase, shortInt, RNG + class(mgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + class(RNG), intent(inout) :: rand + integer(shortInt) :: G + end function sampleEnergyGroup + + !! + !! Sample the time taken for a material particle to transform into a photon + !! Used for ISMC only + !! + function sampleTransformTime(self, matIdx, rand) result(t) + import :: mgIMCDatabase, shortInt, RNG, defReal + class(mgIMCDatabase), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + class(RNG), intent(inout) :: rand + real(defReal) :: t + end function sampleTransformTime + + 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..b1166d9ee --- /dev/null +++ b/NuclearData/mgIMCData/mgIMCMaterial_inter.f90 @@ -0,0 +1,167 @@ +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(getFleck), deferred :: getFleck + procedure(getTemp), deferred :: getTemp + + 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 + + !! + !! 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 + + 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 8b16e96e7..546d34329 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, P_MATERIAL_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 @@ -72,6 +74,9 @@ module nuclearDataReg_mod ! Neutron MG use baseMgNeutronDatabase_class, only : baseMgNeutronDatabase + ! Photon MG + use baseMgIMCDatabase_class, only : baseMgIMCDatabase + implicit none private @@ -97,6 +102,7 @@ module nuclearDataReg_mod public :: kill public :: getNeutronCE public :: getNeutronMG + public :: getIMCMG public :: get public :: getMatNames @@ -110,7 +116,8 @@ module nuclearDataReg_mod !! Parameters character(nameLen), dimension(*), parameter :: AVAILABLE_NUCLEAR_DATABASES = & ['aceNeutronDatabase ', & - 'baseMgNeutronDatabase '] + 'baseMgNeutronDatabase ', & + 'baseMgIMCDatabase '] !! Members type(ndBox),dimension(:),allocatable,target :: databases @@ -122,6 +129,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 !! @@ -318,6 +328,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)) @@ -353,11 +370,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) @@ -393,6 +417,10 @@ subroutine kill() activeIdx_mgNeutron = 0 active_mgNeutron => null() + ! MG IMC + activeIdx_mgIMC = 0 + active_mgIMC => null() + end subroutine kill !! @@ -415,13 +443,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 @@ -433,6 +461,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 !! @@ -459,6 +506,13 @@ function get_byType(type, where) result(ptr) case(P_NEUTRON_MG) ptr => getNeutronMG() + case(P_PHOTON_MG) + ptr => getIMCMG() + + case(P_MATERIAL_MG) + ! Currently only used for ISMC so point to same database + ptr => getIMCMG() + case default ptr => null() end select @@ -572,6 +626,9 @@ subroutine new_nuclearDatabase(database, type) case('baseMgNeutronDatabase') allocate(baseMgNeutronDatabase :: database) + case('baseMgIMCDatabase') + allocate(baseMgIMCDatabase :: database) + case default ! Print available nuclear database types print '(A)', "<><><><><><><><><><><><><><><><><><><><>" 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 17a73d167..0fd52efa7 100644 --- a/ParticleObjects/Source/CMakeLists.txt +++ b/ParticleObjects/Source/CMakeLists.txt @@ -1,8 +1,10 @@ # 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 materialSource_class.f90 + imcMaterialSource_class.f90 + blackBodySource_class.f90 ) diff --git a/ParticleObjects/Source/blackBodySource_class.f90 b/ParticleObjects/Source/blackBodySource_class.f90 new file mode 100644 index 000000000..4d24790b0 --- /dev/null +++ b/ParticleObjects/Source/blackBodySource_class.f90 @@ -0,0 +1,525 @@ +module blackBodySource_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 + use simulationTime_class + + use energyGrid_class, only : energyGrid + use energyGridRegistry_mod, only : get_energyGrid + + implicit none + private + + ! Options for source distribution + integer(shortInt), parameter, public :: SURFACE = 1 + integer(shortInt), parameter, public :: UNIFORM = 2 + integer(shortInt), parameter, public :: OLSON1D = 3 + + !! + !! Generates a source representing a black body + !! Standard is a surface distribution but can be configured for custom distributions + !! + !! 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 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 + !! sampleWeight -> set particle energy-weight + !! sampleTime -> set particle time + !! kill -> terminate source + !! + !! Sample Dictionary Input: + !! source { + !! type blackBodySource; + !! distribution surface; + !! surface -x; + !! temp 1; -> temperature of the black body source + !! } + !! + !! Current source distributions: + !! surface -> black body surface source placed on the surface indicated + !! uniform -> uniform in space, isotropic + !! olson1D -> 1D multifrequency benchmark source from "Stretched and Filtered Multigroup Pn + !! Transport for Improved Positivity and Accuracy", Olson, Gordon Lee, 2020 + !! See materialEquations.f90 for sigma and cv equations for this benchmark. + !! + type, public,extends(configSource) :: blackBodySource + private + + ! Settings defining sourced position and direction for surface sources + real(defReal), dimension(3) :: r = ZERO ! Corner position of source + real(defReal), dimension(3) :: dr = ZERO ! Spatial extent from corner + integer(shortInt), dimension(3,3) :: rotation = ZERO ! Direction rotation matrix + ! Other settings + integer(shortInt) :: distribution = SURFACE + integer(shortInt) :: particleType = P_PHOTON + logical(defBool) :: isMG = .true. + real(defReal) :: T = ZERO ! Source temperature + real(defReal) :: particleWeight = ZERO ! Weight of each particle (= sourceWeight/N) + integer(shortInt) :: timeStepMax = 0 ! Time step to switch source off + ! Probability of emission from each energy group for multi-frequency case + real(defReal), dimension(:), allocatable :: cumEnergyProbs + + contains + ! Initialisation procedures + procedure :: init + procedure :: initSurface + procedure :: initCustom + ! Sampling procedures + procedure :: append + procedure :: sampleType + procedure :: samplePosition + procedure :: sampleEnergy + procedure :: sampleEnergyAngle + procedure :: sampleWeight + procedure :: sampleTime + procedure :: kill + + end type blackBodySource + +contains + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Particle sampling procedures +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! 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(blackBodySource), 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 (blackBodySource_class.f90)' + + ! Store N for calculation of each particle weight + self % particleWeight = self % sourceWeight / N + + ! Generate N particles to populate dungeon + !$omp parallel + pRand = rand + !$omp do private(pRand) + do i = 1, N + call pRand % stride(i) + call dungeon % detain(self % sampleParticle(pRand)) + end do + !$omp end do + !$omp end parallel + + ! Turn off for next time step if needed + if (thisStep() == self % timeStepMax) self % sourceWeight = ZERO + + end subroutine append + + !! + !! Provide particle type + !! + !! See configSource_inter for details. + !! + subroutine sampleType(self, p, rand) + class(blackBodySource), 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(blackBodySource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + integer(shortInt) :: i + real(defReal), dimension(3) :: r + real(defReal) :: x + character(100), parameter :: Here = 'samplePosition (blackBodySource_class.f90)' + + select case(self % distribution) + + case(SURFACE, UNIFORM) + ! Set new x, y and z coords + do i = 1, 3 + r(i) = self % r(i) + rand % get()*self % dr(i) + end do + + case(OLSON1D) + ! Q(x) proportional to exp(-693x**3) (integral from 0 to 4.8 = 0.100909) + rejection:do + x = rand % get() * 4.8_defReal + if (rand % get() < exp(-693*x**3)/0.100909) exit + end do rejection + r(1) = x - 2.4_defReal + r(2) = rand % get() - 0.5_defReal + r(3) = rand % get() - 0.5_defReal + + case default + call fatalError(Here, 'Unrecognised source distribution') + + end select + + ! Assign to particle + p % r = r + + end subroutine samplePosition + + !! + !! Sample angle + !! + !! See configSource_inter for details. + !! + subroutine sampleEnergyAngle(self, p, rand) + class(blackBodySource), 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 (blackBodySource_class.f90)' + + ! Sample direction for isotropic source + if (all(self % rotation == ZERO)) then + ! Sample uniformly within 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 to particle + p % dir = dir + return + end if + + ! If not isotropic (e.g. surface distribution), sample first with a primary direction of +x + phi = TWO_PI * rand % get() + mu = sqrt(rand % get()) + dir = [mu, sqrt(1-mu*mu)*cos(phi), sqrt(1-mu*mu)*sin(phi)] + + ! Rotate to direction of surface normal + p % dir = matmul(self % rotation, dir) + + end subroutine sampleEnergyAngle + + !! + !! Provide particle energy + !! + !! See configSource_inter for details. + !! + subroutine sampleEnergy(self, p, rand) + class(blackBodySource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + real(defReal) :: random + integer(shortInt) :: i + character(100), parameter :: Here = 'sampleEnergy (blackBodySource_class.f90)' + + p % isMG = .true. + + ! Grey case + if (.not.allocated(self % cumEnergyProbs)) then + p % G = 1 + return + end if + + ! Sample energy group + random = rand % get() + do i=1, size(self % cumEnergyProbs) + if (random <= self % cumEnergyProbs(i)) then + p % G = i + return + end if + end do + + call fatalError(Here, 'Somehow failed to sample particle energy group') + + end subroutine sampleEnergy + + !! + !! Provide particle energy-weight + !! + !! See configSource_inter for details. + !! + subroutine sampleWeight(self, p, rand) + class(blackBodySource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + p % wgt = self % particleWeight + + end subroutine sampleWeight + + !! + !! Sample time uniformly within time step + !! + subroutine sampleTime(self, p, rand) + class(blackBodySource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + ! Sample time uniformly within time step + p % time = time % stepStart + timeStep() * rand % get() + + end subroutine sampleTime + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Source initialisation procedures +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Initialise from dictionary + !! + !! See source_inter for details + !! + subroutine init(self, dict, geom) + class(blackBodySource), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), pointer, intent(in) :: geom + integer(shortInt) :: i, j, nGroups + real(defReal) :: nu, eStep + type(energyGrid) :: eGrid + logical(defBool) :: err + character(nameLen) :: gridName, distribution + character(100), parameter :: Here = 'init (blackBodySource_class.f90)' + + ! Provide geometry info to source + self % geom => geom + + ! Provide particle type + self % particleType = P_PHOTON + + ! Get source temperature + call dict % get(self % T, 'temp') + + ! Initialise specifics of source (e.g. position samlping bounds) + call dict % getOrDefault(distribution, 'distribution', 'surface') + select case(distribution) + case('surface') + call self % initSurface(dict) + case default + call self % initCustom(dict) + end select + + ! Get time step to turn off source + call dict % getOrDefault(self % timeStepMax, 'untilStep', 0) + + ! Exit in grey case + gridName = 'mg' + call get_energyGrid(eGrid, gridName, err) + if (err .eqv. .true.) return + + ! Calculate emission probability in each energy group + nGroups = eGrid % getSize() + allocate(self % cumEnergyProbs(nGroups)) + + do i=1, nGroups + eStep = (eGrid % bin(i) - eGrid % bin(i+1)) / 1000 + nu = eGrid % bin(i) - eStep/2 + ! Add previous group probability to get cumulative distribution + if (i > 1) self % cumEnergyProbs(i) = self % cumEnergyProbs(i-1) + ! Step through energies + do j=1, 1000 + self % cumEnergyProbs(i) = self % cumEnergyProbs(i) + eStep*normPlanck(nu,self % T) + nu = nu - eStep + end do + end do + + ! Normalise to account for exclusion of tail ends of distribution + ! TODO: should possibly add these tails into outer energy groups rather than normalising over all groups? + self % cumEnergyProbs = self % cumEnergyProbs / self % cumEnergyProbs(nGroups) + + end subroutine init + + + !! + !! Initialise source for standard black body surface by placing source as one side of + !! bounding box of geometry + !! + !! Input dict should contain 'surface', corresponding to which side of the box is the source + !! e.g. surface -x; => source placed on negative x side of bounding box + !! surface +z; => source placed on positive z side of bounding box + !! etc. + !! + !! Automatically sets bounds for position sampling, rotation matrix for direction sampling, and + !! total weight emitted from the source during each time step + !! + subroutine initSurface(self, dict) + class(blackBodySource), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal), dimension(6) :: bounds + character(nameLen) :: whichSurface + integer(shortInt) :: i + integer(shortInt), dimension(9) :: rotation + real(defReal) :: area + character(100), parameter :: Here = 'initSurface (blackBodySource_class.f90)' + + self % distribution = SURFACE + + bounds = self % geom % bounds() + + ! Bottom left corner + self % r = bounds(1:3) + ! Dimensions of bounding box + do i = 1, 3 + self % dr(i) = bounds(i+3) - bounds(i) + end do + + ! Get position bounds + call dict % get(whichSurface, 'surface') + + select case(whichSurface) + case('-x') + ! Set sampling position to be at constant x value + self % dr(1) = ZERO + ! Nudge to ensure sourcing in correct material + self % r(1) = bounds(1) + TWO*SURF_TOL + ! Set rotation matrix for direction sampling + rotation = [[1,0,0],[0,1,0],[0,0,1]] + case('-y') + self % dr(2) = ZERO + self % r(2) = bounds(2) + TWO*SURF_TOL + rotation = [[0,1,0],[1,0,0],[0,0,1]] + case('-z') + self % dr(3) = ZERO + self % r(3) = bounds(3) + TWO*SURF_TOL + rotation = [[0,0,1],[0,1,0],[1,0,0]] + case('+x') + self % r(1) = bounds(4) - TWO*SURF_TOL + self % dr(1) = ZERO + rotation = [[-1,0,0],[0,1,0],[0,0,1]] + case('+y') + self % r(2) = bounds(5) - TWO*SURF_TOL + self % dr(2) = ZERO + rotation = [[0,-1,0],[1,0,0],[0,0,1]] + case('+z') + self % r(3) = bounds(6) - TWO*SURF_TOL + self % dr(3) = ZERO + rotation = [[0,0,-1],[0,1,0],[1,0,0]] + case default + call fatalError(Here, 'Unrecognised surface, must be +/-x,y or z') + end select + + ! Note that reshape fills columns first so the above rotations are [[col1][col2][col3]] + self % rotation = reshape(rotation,[3,3]) + + ! Calculate surface area of source + area = product(self % dr, self % dr /= ZERO) + + ! Calculate total source energy + self % sourceWeight = radiationConstant * lightSpeed * timeStep() * self % T**4 * area / 4 + + end subroutine initSurface + + !! + !! Initialise black body source for more unique cases, e.g. Olson 1D MG benchmark (2020) + !! + !! Rotation matrix of ZERO corresponds to isotropic direction sampling. Set source weight or + !! other settings as required for specific case, and add additional cases to individual samping + !! procedures if needed. + !! + subroutine initCustom(self, dict) + class(blackBodySource), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(nameLen) :: name + integer(shortInt) :: i + real(defReal), dimension(6) :: bounds + character(100), parameter :: Here = 'initCustom (blackBodySource_class.f90)' + + call dict % get(name, 'distribution') + + select case(name) + + case('uniform') + self % distribution = UNIFORM + bounds = self % geom % bounds() + ! Bottom left corner + self % r = bounds(1:3) + ! Dimensions of bounding box + do i = 1, 3 + self % dr(i) = bounds(i+3) - bounds(i) + end do + ! Isotropic direction sampling + self % rotation = ZERO + call dict % get(self % sourceWeight, 'sourceWeight') + + case('olson1D') + ! See self % samplePosition for position sampling specifics for Olson1D + self % distribution = OLSON1D + ! Isotropic directional sampling + self % rotation = ZERO + ! Set source weight + self % sourceWeight = radiationConstant * lightSpeed * timeStep() * self % T**4 * 0.100909 + + case default + call fatalError(Here, 'Unrecognised distribution for custom black body source') + + end select + + end subroutine initCustom + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(blackBodySource), intent(inout) :: self + + ! Kill superclass + call kill_super(self) + + ! Kill local components + self % r = ZERO + self % dr = ZERO + self % distribution = SURFACE + self % particleType = P_PHOTON + self % isMG = .true. + self % T = ZERO + self % particleWeight = ZERO + + end subroutine kill + + +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!! Misc. procedures +!!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + + !! + !! Return the frequency-normalised Planck spectrum evaluated for given frequency + !! nu and temperature T + !! + pure function normPlanck(E, T) result(b) + real(defReal), intent(in) :: E, T + real(defReal) :: b + + b = 15*E**3 / ((pi*T)**4 * (exp(E/T)-1)) + + end function normPlanck + + +end module blackBodySource_class diff --git a/ParticleObjects/Source/configSource_inter.f90 b/ParticleObjects/Source/configSource_inter.f90 index 968dd0dbe..0436b2530 100644 --- a/ParticleObjects/Source/configSource_inter.f90 +++ b/ParticleObjects/Source/configSource_inter.f90 @@ -32,6 +32,8 @@ module configSource_inter contains procedure :: sampleParticle + procedure :: sampleWeight + procedure :: sampleTime procedure(sampleType), deferred :: sampleType procedure(samplePosition), deferred :: samplePosition procedure(sampleEnergy), deferred :: sampleEnergy @@ -133,11 +135,37 @@ function sampleParticle(self, rand) result(p) call self % samplePosition(p, rand) call self % sampleEnergyAngle(p, rand) call self % sampleEnergy(p, rand) - p % time = ZERO - p % wgt = ONE + call self % sampleWeight(p, rand) + call self % sampleTime(p, rand) 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 + + !! + !! Set particle's time to 0 + !! Can be overriden in subclasses if needed + !! + subroutine sampleTime(self, p, rand) + class(configSource), intent(inout) :: self + class(particleState), intent(inout) :: p + class(RNG), intent(inout) :: rand + + p % time = ZERO + + end subroutine sampleTime + !! !! Return to uninitialised state !! diff --git a/ParticleObjects/Source/imcMaterialSource_class.f90 b/ParticleObjects/Source/imcMaterialSource_class.f90 new file mode 100644 index 000000000..ef8cf4e17 --- /dev/null +++ b/ParticleObjects/Source/imcMaterialSource_class.f90 @@ -0,0 +1,284 @@ +module imcMaterialSource_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, P_MATERIAL + use particleDungeon_class, only : particleDungeon + use source_inter, only : source, kill_super => kill + + use geometry_inter, only : geometry + use geometryGrid_class, only : geometryGrid + 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 + + use simulationTime_class + + implicit none + private + + ! Calculation type + integer(shortInt), parameter :: IMC = 1, ISMC = 2 + + !! + !! Source for uniform generation of photons within a material for IMC and ISMC calculations + !! + !! 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 -> Energy group + !! pType -> P_PHOTON for IMC, P_MATERIAL for ISMC + !! bounds -> Bounds of geometry + !! calcType -> IMC or ISMC, changes type of material to be sampled + !! + !! Interface: + !! source_inter Interface + !! + !! SAMPLE INPUT: + !! matSource { type imcMaterialSource; calcType IMC; } + !! + type, public,extends(source) :: imcMaterialSource + private + logical(defBool) :: isMG = .true. + real(defReal), dimension(3) :: bottom = ZERO + real(defReal), dimension(3) :: top = ZERO + integer(shortInt) :: G = 0 + integer(shortInt) :: pType = P_PHOTON + real(defReal), dimension(6) :: bounds = ZERO + integer(shortInt) :: calcType = IMC + contains + procedure :: init + procedure :: append + procedure :: sampleParticle + procedure, private :: sampleIMC + procedure :: kill + end type imcMaterialSource + +contains + + !! + !! Initialise material Source + !! + !! See source_inter for details + !! + subroutine init(self, dict, geom) + class(imcMaterialSource), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), pointer, intent(in) :: geom + character(100), parameter :: Here = 'init (imcMaterialSource_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 calculation type - Automatically added to dict in implicitPhysicsPackage + call dict % getOrDefault(self % calcType, 'calcType', IMC) + select case(self % calcType) + case(IMC) + self % pType = P_PHOTON + case(ISMC) + self % pType = P_MATERIAL + case default + call fatalError(Here, 'Unrecognised calculation type. Should be "IMC" or "ISMC"') + end select + + end subroutine init + + + !! + !! Generate N particles to add to a particleDungeon without overriding particles already present. + !! Note that energy here refers to energy weight rather than group. + !! + !! 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(imcMaterialSource), 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, G + real(defReal) :: energy, totalEnergy + type(RNG) :: pRand + class(mgIMCDatabase), pointer :: nucData + class(geometry), pointer :: geom + character(100), parameter :: Here = "append (imcMaterialSource_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 + if (self % calcType == IMC) then + totalEnergy = nucData % getEmittedRad() + else + totalEnergy = nucData % getMaterialEnergy() + end if + + ! Loop through materials + do matIdx = 1, mm_nMat() + + ! Get energy to be emitted from material matIdx + if (self % calcType == IMC) then + energy = nucData % getEmittedRad(matIdx) + else + energy = nucData % getMaterialEnergy(matIdx) + end if + + ! 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 + geom => self % geom + select type(geom) + class is(geometryGrid) + bounds = geom % matBounds(matIdx) + class default + bounds = self % bounds + end select + + ! Find energy per particle + energy = energy / Ntemp + + ! Sample particles + !$omp parallel + pRand = rand + !$omp do private(pRand, G) + do i=1, Ntemp + call pRand % stride(i) + G = nucData % sampleEnergyGroup(matIdx, pRand) + call dungeon % detain(self % sampleIMC(pRand, matIdx, energy, G, bounds)) + end do + !$omp end do + !$omp end parallel + + end if + end do + + end subroutine append + + !! + !! Should not be called + !! + function sampleParticle(self, rand) result(p) + class(imcMaterialSource), intent(inout) :: self + class(RNG), intent(inout) :: rand + type(particleState) :: p + character(100), parameter :: Here = 'sampleParticle (imcMaterialSource_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 + !! G [in] -> energy group of sampled particle + !! bounds [in] -> bounds for position search, will be bounds of entire geometry if using + !! geometryStd, and bounds of single material if using geometryGrid + !! + function sampleIMC(self, rand, targetMatIdx, energy, G, bounds) result(p) + class(imcMaterialSource), intent(inout) :: self + class(RNG), intent(inout) :: rand + integer(shortInt), intent(in) :: targetMatIdx + real(defReal), intent(in) :: energy + integer(shortInt), intent(in) :: G + 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 (imcMaterialSource_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 + + ! 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) + + ! Sample time uniformly within time step + p % time = time % stepStart + timeStep() * rand % get() + + ! Assign basic phase-space coordinates + p % matIdx = matIdx + p % uniqueID = uniqueID + p % r = r + p % dir = dir + p % G = G + p % isMG = .true. + p % wgt = energy + p % type = self % pType + + end function sampleIMC + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(imcMaterialSource), intent(inout) :: self + + call kill_super(self) + + self % isMG = .true. + self % bounds = ZERO + self % G = 0 + + end subroutine kill + +end module imcMaterialSource_class diff --git a/ParticleObjects/Source/sourceFactory_func.f90 b/ParticleObjects/Source/sourceFactory_func.f90 index 474001979..c279a694d 100644 --- a/ParticleObjects/Source/sourceFactory_func.f90 +++ b/ParticleObjects/Source/sourceFactory_func.f90 @@ -8,9 +8,11 @@ module sourceFactory_func use source_inter, only : source ! source implementations - use pointSource_class, only : pointSource - use fissionSource_class, only : fissionSource - use materialSource_class, only : materialSource + use pointSource_class, only : pointSource + use fissionSource_class, only : fissionSource + use materialSource_class, only : materialSource + use blackBodySource_class, only : blackBodySource + use imcMaterialSource_class, only : imcMaterialSource ! geometry use geometry_inter, only : geometry @@ -24,9 +26,11 @@ 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 ',& - 'materialSource'] + character(nameLen),dimension(*),parameter :: AVAILABLE_sources = [ 'pointSource ',& + 'fissionSource ',& + 'materialSource ',& + 'blackBodySource ',& + 'imcMaterialSource'] contains @@ -58,6 +62,12 @@ subroutine new_source(new, dict, geom) case('materialSource') allocate(materialSource :: new) + case('blackBodySource') + allocate(blackBodySource :: new) + + case('imcMaterialSource') + allocate(imcMaterialSource :: new) + case default print *, AVAILABLE_sources call fatalError(Here, 'Unrecognised type of source: ' // trim(type)) diff --git a/ParticleObjects/Source/source_inter.f90 b/ParticleObjects/Source/source_inter.f90 index 141cc3558..66ab3e66c 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,14 +30,17 @@ 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 !! type, public,abstract :: source private class(geometry), pointer, public :: geom => null() + real(defReal), public :: sourceWeight = ZERO contains - procedure, non_overridable :: generate + procedure :: generate + procedure :: append procedure(sampleParticle), deferred :: sampleParticle procedure(init), deferred :: init procedure(kill), deferred :: kill @@ -122,6 +126,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 7d47337ec..784d9ac0c 100644 --- a/ParticleObjects/particleDungeon_class.f90 +++ b/ParticleObjects/particleDungeon_class.f90 @@ -1,9 +1,11 @@ module particleDungeon_class use numPrecision - use genericProcedures, only : fatalError, numToChar - use particle_class, only : particle, particleState + use genericProcedures, only : fatalError, numToChar, linFind, getDistance + use particle_class, only : particle, particleState, P_MATERIAL, P_PHOTON use RNG_class, only : RNG + use geometry_inter, only : geometry + use universalVariables, only : INF implicit none private @@ -46,22 +48,28 @@ module particleDungeon_class !! get(i) -> function returns particle state at index i !! !! Misc procedures: - !! isEmpty() -> returns .true. if there are no more particles - !! cleanPop() -> kill or prisoners - !! normWeight(totWgt)-> normalise dungeon population so its total weight is totWgt - !! normSize(N) -> normalise dungeon population so it contains N particles - !! 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" + !! isEmpty() -> returns .true. if there are no more particles + !! cleanPop() -> kill or prisoners + !! normWeight(totWgt) -> normalise dungeon population so its total weight is totWgt + !! normSize(N) -> normalise dungeon population so it contains N particles + !! does not take nonuniform weight of particles into account + !! combine(idx1,idx2) -> combine 2 particles by summing their weight and moving to a weight- + !! averaged position + !! deleteParticle(idx) -> deletes particle at idx and reduces dungeon size by 1 + !! 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 - !! kill() -> return to uninitialised state + !! init(maxSize) -> allocate space to store maximum of maxSize particles + !! kill() -> return to uninitialised state !! 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, public :: prisoners @@ -86,11 +94,15 @@ module particleDungeon_class procedure :: isEmpty procedure :: normWeight procedure :: normSize + procedure :: closest + procedure :: combine + procedure :: deleteParticle procedure :: cleanPop procedure :: popSize procedure :: popWeight procedure :: setSize procedure :: printToFile + procedure :: printToScreen ! Private procedures procedure, private :: detain_particle @@ -402,7 +414,7 @@ subroutine normSize(self,N,rand) class(RNG), intent(inout) :: rand integer(shortInt) :: excessP integer(shortInt) :: i, idx - character(100), parameter :: Here =' normSize (particleDungeon_class.f90)' + character(100), parameter :: Here = 'normSize (particleDungeon_class.f90)' ! Protect against invalid N if( N > size(self % prisoners)) then @@ -436,6 +448,108 @@ subroutine normSize(self,N,rand) end subroutine normSize + !! + !! Find the closest particle to particle at idx, that is of the same type and in the same material + !! + function closest(self, idx) result(idxClose) + class(particleDungeon), intent(in) :: self + integer(shortInt), intent(in) :: idx + integer(shortInt) :: idxClose, matIdx, type, i + real(defReal), dimension(3) :: r + real(defReal) :: dist, minDist + character(100), parameter :: Here = 'closest (particleDungeon_class.f90)' + + ! Get required properties of particle + r = self % prisoners(idx) % r + type = self % prisoners(idx) % type + matIdx = self % prisoners(idx) % matIdx + + minDist = INF + idxClose = 0 + + !$omp parallel + !$omp do private(dist) + do i=1, self % pop + ! Require particles to be of same type and in same matIdx + if (self % prisoners(i) % matIdx /= matIdx) cycle + if (self % prisoners(i) % type /= type) cycle + + ! Get distance + dist = getDistance(r, self % prisoners(i) % r) + if (dist < minDist) then + minDist = dist + idxClose = i + end if + + end do + !$omp end do + !$omp end parallel + + if (idxClose == 0) call fatalError(Here, 'No valid particle found') + + end function closest + + !! + !! Combine two particles in the dungeon by summing their weight and moving to a weighted- + !! average position. Direction is unchanged. + !! + !! Particle at idx1 is moved to a position that is the energy-weighted average + !! of the two original positions. Its new energy is the sum of the two original energies. + !! + !! Particle at idx2 remains unchanged. To delete, call self % deleteParticle(idx2) afterwards. + !! + !! Args: + !! idx1 [in] -> Index of 1st particle, will be overridden + !! idx2 [in] -> Index of 2nd particle, will be kept unchanged + !! + subroutine combine(self, idx1, idx2) + class(particleDungeon), intent(inout) :: self + integer(shortInt), intent(in) :: idx1 + integer(shortInt), intent(in) :: idx2 + type(particle) :: p1, p2 + real(defReal), dimension(3) :: r1, r2, rNew + + ! Get initial particle data + call self % copy(p1, idx1) + call self % copy(p2, idx2) + r1 = p1 % rGlobal() + r2 = p2 % rGlobal() + + ! Move to new combined position + rNew = (r1*p1 % w + r2*p2 % w) / (p1 % w + p2 % w) + call p1 % teleport(rNew) + + ! Combine weights and overwrite particle + p1 % w = p1 % w + p2 % w + call self % replace(p1, idx1) + + end subroutine combine + + !! + !! Deletes particle at idx by releasing top particle and copying into position idx, + !! overwriting particle to be deleted + !! + !! Args: + !! idx [in] -> index of particle to be deleted + !! + subroutine deleteParticle(self, idx) + class(particleDungeon), intent(inout) :: self + integer(shortInt), intent(in) :: idx + type(particle) :: p + integer(shortInt) :: matIdx + + matIdx = self % prisoners(self % pop) % matIdx + + ! Release particle at top of dungeon + call self % release(p) + + ! Copy into position of particle to be deleted + if (idx /= self % pop + 1) call self % replace(p, idx) + + self % prisoners(idx) % matIdx = matIdx + + end subroutine deleteParticle + !! !! Kill or particles in the dungeon !! @@ -447,7 +561,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 +631,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 +647,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 0fde0fd91..a3c10dc98 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -1,5 +1,6 @@ module particle_class + use numPrecision use universalVariables use genericProcedures @@ -12,8 +13,9 @@ module particle_class !! !! Particle types paramethers !! - integer(shortInt), parameter,public :: P_NEUTRON = 1,& - P_PHOTON = 2 + integer(shortInt), parameter,public :: P_NEUTRON = 1,& + P_PHOTON = 2,& + P_MATERIAL = 3 !! !! Public particle type procedures @@ -138,6 +140,7 @@ module particle_class procedure :: getUniIdx procedure :: matIdx procedure, non_overridable :: getType + procedure :: getSpeed ! Operations on coordinates procedure :: moveGlobal @@ -402,13 +405,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, P_MATERIAL_MG !! !! Errors: !! None @@ -417,14 +420,51 @@ 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 + + ! Currently only MG material particles supported + else if (self % type == P_MATERIAL) then + type = P_MATERIAL_MG 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 !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> @@ -709,7 +749,6 @@ elemental subroutine kill_particleState(self) end subroutine kill_particleState - !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> !! Misc Procedures !!<><><><><><><>><><><><><><><><><><><>><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> @@ -726,6 +765,7 @@ elemental function verifyType(type) result(isValid) ! Check against particles types isValid = isValid .or. type == P_NEUTRON isValid = isValid .or. type == P_PHOTON + isValid = isValid .or. type == P_MATERIAL end function verifyType @@ -743,6 +783,9 @@ pure function printType(type) result(name) case(P_PHOTON) name = 'Photon' + case(P_MATERIAL) + name = 'Material' + case default name = 'INVALID PARTICLE TYPE' diff --git a/PhysicsPackages/CMakeLists.txt b/PhysicsPackages/CMakeLists.txt index a1db7a1f4..8ddb23af6 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 + ./implicitPhysicsPackage_class.f90 ./vizPhysicsPackage_class.f90 ./rayVolPhysicsPackage_class.f90 ) diff --git a/PhysicsPackages/implicitPhysicsPackage_class.f90 b/PhysicsPackages/implicitPhysicsPackage_class.f90 new file mode 100644 index 000000000..f83092b8c --- /dev/null +++ b/PhysicsPackages/implicitPhysicsPackage_class.f90 @@ -0,0 +1,597 @@ +module implicitPhysicsPackage_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, P_MATERIAL + 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_geomIdx => geomIdx + use geometryFactory_func, only : new_geometry + 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 energyWeightClerk_class, only : energyWeightClerkResult + + ! Factories + use transportOperatorFactory_func, only : new_transportOperator + use sourceFactory_func, only : new_source + + use simulationTime_class + use energyGridRegistry_mod, only : define_energyGrid + + implicit none + + private + + ! Calculation Type + integer(shortInt), parameter, public :: IMC = 1 + integer(shortInt), parameter, public :: ISMC = 2 + + !! + !! Physics Package for Implicit Monte Carlo calculations + !! + !! Settings: + !! method -> IMC or ISMC + !! pop -> For IMC, approx. number of new particles to generate per time step + !! -> For ISMC, starting population of material particles + !! limit -> Size of particle dungeons + !! steps -> Number of time steps to simulate + !! timeStep -> Time step to be used + !! printUpdates (OPTIONAL) -> Prints material update info for first N materials + !! + type, public,extends(physicsPackage) :: implicitPhysicsPackage + private + ! Building blocks + 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 :: energyWeightAtch => null() + + ! Settings + integer(shortInt) :: N_steps + integer(shortInt) :: pop + integer(shortInt) :: limit + integer(shortInt) :: method + character(pathLen) :: outputFile + character(nameLen) :: outputFormat + integer(shortInt) :: printSource = 0 + 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 :: matSource + + ! 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 implicitPhysicsPackage + +contains + + subroutine run(self) + class(implicitPhysicsPackage), intent(inout) :: self + + print *, repeat("<>",50) + print *, "/\/\ IMPLICIT CALCULATION /\/\" + + call self % steps(self % tally, self % energyWeightAtch, self % N_steps) + call self % collectResults() + + print * + print *, "\/\/ END OF IMPLICIT CALCULATION \/\/" + print * + end subroutine + + !! + !! Run steps for calculation + !! + !! + !! Notes differences between IMC and ISMC regarding particle generation: + !! + !! -> IMC generates particles from material emission as well as input source, ISMC only + !! generates from input source. + !! + !! -> Particles for IMC are killed when absorbed, but for ISMC remain as material particles. + !! This allows IMC to keep dungeon pop under control far better than ISMC. For ISMC we + !! generate particles such that pop is at limit at end of calculation, with runtime per + !! time step increasing approximatelty linearly as the calculation progresses + !! + subroutine steps(self, tally, tallyAtch, N_steps) + class(implicitPhysicsPackage), 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, nFromMat, num, nParticles + type(particle), save :: p + real(defReal) :: sourceWeight, elapsed_T, end_T, T_toEnd + class(IMCMaterial), pointer :: mat + character(100),parameter :: Here ='steps (implicitPhysicsPackage_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) + + ! Generate starting population of material particles for ISMC + if (self % method == ISMC) then + call self % matSource % append(self % thisStep, self % pop, self % pRNG) + nFromMat = 0 + end if + + do i=1,N_steps + + ! Generate particles while staying below dungeon limit (see note in subroutine description) + if (self % method == IMC) then + + ! Reduce number of particles to generate if close to limit + 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 + N = max(1, N) + end if + + ! Calculate proportion to be generated from input source + if (allocated(self % inputSource)) then + sourceWeight = self % inputSource % sourceWeight + nFromMat = int(N * (1 - sourceWeight/(sourceWeight + self % nucData % getEmittedRad()))) + ! Generate from input source + call self % inputSource % append(self % thisStep, N - nFromMat, self % pRNG) + end if + + ! Add to dungeon particles emitted from material + call self % matSource % append(self % thisStep, nFromMat, self % pRNG) + + ! ISMC particle generation + else if (allocated(self % inputSource)) then + + ! Generate particles such that pop is almost at limit at calculation end + N = (self % limit - self % thisStep % popSize()) / (N_steps-i+1) + call self % inputSource % append(self % thisStep, N, 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 .and. p % getType() /= P_MATERIAL_MG) then + call fatalError(Here, 'Particle is not of type P_PHOTON_MG or P_MATERIAL_MG') + end if + + ! Assign maximum particle time + p % timeMax = time % stepEnd + + ! Check for time errors + if (p % time < time % stepStart .or. p % time >= time % stepEnd) then + call fatalError(Here, 'Particle time not within time step') + 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) + + ! Cycle if particle history not yet completed + if (self % method == ISMC .or. p % type == P_PHOTON) cycle history + + ! If P_MATERIAL and IMC, kill particle and exit + p % isDead = .true. + p % fate = ABS_FATE + call tally % reportHist(p) + 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, 'energyWeightTally') + + ! Update material properties using tallied energy + select type(tallyRes) + class is(energyWeightClerkResult) + call self % nucData % updateProperties(tallyRes % materialEnergy, self % printUpdates) + class default + call fatalError(Here, 'Tally result class should be energyWeightClerkResult') + end select + + ! Reset tally for next time step + if (i /= N_Steps) call tallyAtch % reset('energyWeightTally') + + ! Advance to next time step + call nextStep() + + ! 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() + + 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) + + ! Output final radiation energies + open(unit = 11, file = 'radEnergy.txt') + select type(tallyRes) + class is(energyWeightClerkResult) + write(11, '(8A)') numToChar(tallyRes % radiationEnergy) + class default + call fatalError(Here, 'Tally result class should be energyWeightClerkResult') + end select + close(11) + + end subroutine steps + + !! + !! Print calculation results to file + !! + subroutine collectResults(self) + class(implicitPhysicsPackage), 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(implicitPhysicsPackage), intent(inout) :: self + class(dictionary), intent(inout) :: dict + class(dictionary), pointer :: tempDict + type(dictionary) :: locDict1, locDict2, locDict3 + 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 + real(defReal) :: timeStep + type(dictionary),target :: newGeom, newData + character(nameLen) :: method, units + character(100), parameter :: Here ='init (implicitPhysicsPackage_class.f90)' + + call cpu_time(self % CPU_time_start) + + ! Get method + call dict % getOrDefault(method, 'method', 'IMC') + select case(method) + case ('IMC') + self % method = IMC + case('ISMC') + self % method = ISMC + case default + call fatalError(Here, 'Unrecognised method') + end select + + ! 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(timeStep,'timeStep') + call dict % getOrDefault(self % printUpdates, 'printUpdates', 0) + nucData = 'mg' + + ! Set time step after changing units if necessary + call dict % getOrDefault(units, 'units', 'ns') + select case(units) + case('s') + ! No change needed + case('ns') + ! Convert time step from ns to s + timeStep = timeStep * 1e-9_defReal + case('marshak') + ! Special case where a = c = 1 + timeStep = timeStep/lightSpeed + print *, 'WARNING: For Marshak wave, still need to manually change radiationConstant to 1' + case default + call fatalError(Here, 'Unrecognised units') + end select + call setStep(timeStep) + + ! 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 caught 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) + + ! Initialise energy grid in multi-frequency case + if (dict % isPresent('energyGrid')) then + tempDict => dict % getDictPtr('energyGrid') + call define_energyGrid(nucData, tempDict) + print *, 'Energy grid defined: ', nucData + end if + + ! Build Nuclear Data + call ndReg_init(dict % getDictPtr('nuclearData')) + + ! Build geometry + geomName = 'IMCGeom' + call new_geometry(dict % getDictPtr('geometry'), geomName) + self % geomIdx = gr_geomIdx(geomName) + self % geom => gr_geomPtr(self % geomIdx) + + ! Activate Nuclear Data *** All materials are active + call ndReg_activate(P_PHOTON_MG, nucData, self % geom % activeMats()) + self % nucData => mgIMCDatabase_CptrCast(ndReg_get(P_PHOTON_MG)) + + call newGeom % kill() + call newData % kill() + + ! Initialise material source + call locDict1 % init(2) + call locDict1 % store('type', 'imcMaterialSource') + ! Tell source if we are using IMC or ISMC + call locDict1 % store('calcType', self % method) + call new_source(self % matSource, locDict1, self % geom) + call locDict1 % kill() + + ! Read external particle source definition + if( dict % isPresent('source') ) then + tempDict => dict % getDictPtr('source') + call new_source(self % inputSource, tempDict, self % geom) + 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 calculation type + call self % nucData % setCalcType(self % method) + + ! 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 energy weight tally attachment + call locDict1 % init(1) + call locDict2 % init(2) + call locDict3 % init(2) + call locDict3 % store('type','materialMap') + call locDict3 % store('materials', [mats]) + call locDict2 % store('type','energyWeightClerk') + call locDict2 % store('map', locDict3) + call locDict1 % store('energyWeightTally', locDict2) + + allocate(self % energyWeightAtch) + call self % energyWeightAtch % init(locDict1) + call self % tally % push(self % energyWeightAtch) + + ! 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(implicitPhysicsPackage), intent(inout) :: self + + ! TODO: This subroutine + + end subroutine kill + + !! + !! Print settings of the physics package + !! + subroutine printSettings(self) + class(implicitPhysicsPackage), 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 implicitPhysicsPackage_class diff --git a/PhysicsPackages/physicsPackageFactory_func.f90 b/PhysicsPackages/physicsPackageFactory_func.f90 index 48de37e1c..8d9d816f5 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 implicitPhysicsPackage_class, only : implicitPhysicsPackage implicit none private @@ -26,6 +26,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',& + 'implicitPhysicsPackage ',& 'vizPhysicsPackage ',& 'rayVolPhysicsPackage '] @@ -49,7 +50,7 @@ function new_physicsPackage(dict) result(new) ! Obtain string that specifies type to be built call dict % get(type,'type') - ! Allocate approperiate subclass of physicsPackage + ! Allocate appropriate subclass of physicsPackage select case(type) case('eigenPhysicsPackage') allocate( eigenPhysicsPackage :: new) @@ -57,6 +58,9 @@ function new_physicsPackage(dict) result(new) case('fixedSourcePhysicsPackage') allocate( fixedSourcePhysicsPackage :: new) + case('implicitPhysicsPackage') + allocate( implicitPhysicsPackage :: new) + case('vizPhysicsPackage') allocate( vizPhysicsPackage :: new) diff --git a/SharedModules/CMakeLists.txt b/SharedModules/CMakeLists.txt index dcd1633ec..8320d6e8d 100644 --- a/SharedModules/CMakeLists.txt +++ b/SharedModules/CMakeLists.txt @@ -12,6 +12,7 @@ add_sources( ./genericProcedures.f90 ./timer_mod.f90 ./charLib_func.f90 ./openmp_func.f90 + ./simulationTime_class.f90 ./errors_mod.f90) add_unit_tests( ./Tests/grid_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 2197c3274..bfa253e0b 100644 --- a/SharedModules/endfConstants.f90 +++ b/SharedModules/endfConstants.f90 @@ -107,6 +107,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/genericProcedures.f90 b/SharedModules/genericProcedures.f90 index 2cbb9d280..e15718d5c 100644 --- a/SharedModules/genericProcedures.f90 +++ b/SharedModules/genericProcedures.f90 @@ -1481,4 +1481,16 @@ subroutine printFishLineR(offset) end subroutine printFishLineR + !! + !! Returns the euclidean distance between two 3D points + !! + function getDistance(r1, r2) result(dist) + real(defReal), dimension(3), intent(in) :: r1 + real(defReal), dimension(3), intent(in) :: r2 + real(defReal) :: dist + + dist = (r2(1) - r1(1))**2 + (r2(2) - r1(2))**2 + (r2(3) - r1(3))**2 + + end function getDistance + end module genericProcedures diff --git a/SharedModules/simulationTime_class.f90 b/SharedModules/simulationTime_class.f90 new file mode 100644 index 000000000..72cd1cfa8 --- /dev/null +++ b/SharedModules/simulationTime_class.f90 @@ -0,0 +1,89 @@ + +module simulationTime_class + + use numPrecision + use genericProcedures, only : fatalError + + implicit none + private + + !! + !! Simple module to keep track of simulation time for time-dependent calculations + !! + !! Allows easy public access to time step and current time, and provides a single interface to + !! change time step for calculations with a variable time step, rather than needing to update the + !! time step separately in all required modules (source, material, etc.) + !! + type, public :: simulationTime + real(defReal) :: step = ONE + real(defReal) :: stepStart = ZERO + real(defReal) :: stepEnd = ONE + integer(shortInt) :: stepsCompleted = 0 + end type simulationTime + + type(simulationTime), public :: time + + public :: setStep + public :: thisStep + public :: nextStep + public :: timeStep + public :: timeLeft + +contains + + !! + !! Set time step + !! + subroutine setStep(dt) + real(defReal), intent(in) :: dt + character(100), parameter :: Here = 'set (timeStep_class.f90)' + + if (dt <= ZERO) call fatalError(Here, 'Time step must be positive') + + time % step = dt + time % stepEnd = dt + + end subroutine setStep + + function thisStep() result(i) + integer(shortInt) :: i + + i = time % stepsCompleted + 1 + + end function thisStep + + !! + !! Advance time by one time step + !! + subroutine nextStep() + + time % stepsCompleted = time % stepsCompleted + 1 + + ! Set step start and end time + time % stepStart = time % step * time % stepsCompleted + time % stepEnd = time % step * (time % stepsCompleted + 1) + + end subroutine nextStep + + !! + !! Return time step size + !! + function timeStep() result(dt) + real(defReal) :: dt + + dt = time % step + + end function timeStep + + !! + !! Return time remaining until end of time step + !! + function timeLeft(t) result(remaining_t) + real(defReal), intent(in) :: t + real(defReal) :: remaining_t + + remaining_t = time % stepEnd - t + + end function timeLeft + +end module simulationTime_class diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 9c133c1ee..d030a1282 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -60,8 +60,11 @@ module universalVariables Z_AXIS = 3 ! Particle Type Enumeration - integer(shortInt), parameter :: P_NEUTRON_CE = 1, & - P_NEUTRON_MG = 2 + integer(shortInt), parameter :: P_NEUTRON_CE = 1, & + P_NEUTRON_MG = 2, & + P_PHOTON_CE = 3, & + P_PHOTON_MG = 4, & + P_MATERIAL_MG = 5 ! Search error codes integer(shortInt), parameter :: valueOutsideArray = -1,& @@ -72,7 +75,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 1b0997906..baf0d0da9 100644 --- a/Tallies/TallyClerks/CMakeLists.txt +++ b/Tallies/TallyClerks/CMakeLists.txt @@ -11,6 +11,7 @@ add_sources(./tallyClerk_inter.f90 ./dancoffBellClerk_class.f90 ./shannonEntropyClerk_class.f90 ./centreOfMassClerk_class.f90 + ./energyWeightClerk_class.f90 ./mgXsClerk_class.f90 ) diff --git a/Tallies/TallyClerks/collisionClerk_class.f90 b/Tallies/TallyClerks/collisionClerk_class.f90 index dd5dd66df..054b5ac4b 100644 --- a/Tallies/TallyClerks/collisionClerk_class.f90 +++ b/Tallies/TallyClerks/collisionClerk_class.f90 @@ -94,7 +94,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/energyWeightClerk_class.f90 b/Tallies/TallyClerks/energyWeightClerk_class.f90 new file mode 100644 index 000000000..6f12fec0e --- /dev/null +++ b/Tallies/TallyClerks/energyWeightClerk_class.f90 @@ -0,0 +1,388 @@ +module energyWeightClerk_class + + use numPrecision + use tallyCodes + use genericProcedures, only : fatalError + use dictionary_class, only : dictionary + use particle_class, only : particle, particleState, P_PHOTON + 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 + + !! + !! Clerk for tallying energy weights. Weight response automatically initialised, does not accept + !! additional responses. + !! + !! Tallies two different energies: + !! materialEnergy (reportHist, weight of particles absorbed into material) + !! radiationEnergy (reportTrans, only for particles with fate=AGED_FATE) + !! + !! 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 DICTIONARY INPUT: + !! + !! myEnergyWeightClerk { + !! type energyWeightClerk; + !! # filter { } # + !! # map { } # + !! } + !! + !! Note that tally % reset('myEnergyWeightClerk') should be called between time steps to avoid + !! energy accumulation across time steps + !! + type, public, extends(tallyClerk) :: energyWeightClerk + private + ! Filter, Map & Vector of Responses + class(tallyFilter), allocatable :: filter + class(tallyMap), allocatable :: map + type(tallyResponseSlot),dimension(:),allocatable :: response + + ! Useful data + integer(shortInt) :: width = 0 + + contains + ! Procedures used during build + procedure :: init + procedure :: kill + procedure :: validReports + procedure :: getSize + + ! File reports and check status -> run-time procedures + procedure :: reportHist + procedure :: reportTrans + + ! Output procedures + procedure :: display + procedure :: print + procedure :: getResult + + end type energyWeightClerk + + !! + !! Result class, gives access to tallied material energy and radiation energy + !! + type, public, extends(tallyResult) :: energyWeightClerkResult + real(defReal), dimension(:), allocatable :: materialEnergy + real(defReal), dimension(:), allocatable :: radiationEnergy + end type energyWeightClerkResult + + !! Relative bin positions of each energy + integer(shortInt), parameter :: MATERIAL_ENERGY = 1 + integer(shortInt), parameter :: RADIATION_ENERGY = 2 + +contains + + !! + !! Initialise clerk from dictionary and name + !! + !! See tallyClerk_inter for details + !! + subroutine init(self, dict, name) + class(energyWeightClerk), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(nameLen), intent(in) :: name + character(nameLen),dimension(:),allocatable :: responseNames + integer(shortInt) :: i + type(dictionary) :: locDict + character(100), parameter :: Here = 'init (energyWeightClerk_class.f90)' + + ! 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 + + ! Call error if responses are given + if( dict % isPresent('response')) then + call fatalError(Here, 'Warning: response not needed for energyWeightClerk') + end if + + ! Initialise weight response automatically + call locDict % init(1) + call locDict % store('type','weightResponse') + allocate(self % response(1)) + call self % response(1) % init(locDict) + + self % width = 2 + + end subroutine init + + !! + !! Return to uninitialised state + !! + elemental subroutine kill(self) + class(energyWeightClerk), 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(energyWeightClerk),intent(in) :: self + integer(shortInt),dimension(:),allocatable :: validCodes + + validCodes = [trans_CODE,hist_CODE] + + end function validReports + + !! + !! Return memory size of the clerk + !! + !! See tallyClerk_inter for details + !! + elemental function getSize(self) result(S) + class(energyWeightClerk), intent(in) :: self + integer(shortInt) :: S + + S = self % width + if(allocated(self % map)) S = S * self % map % bins(0) + + end function getSize + + !! + !! Tally an increase in material energy after a particle has been absorbed + !! + subroutine reportHist(self, p, xsData, mem) + class(energyWeightClerk), 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 (energyWeightClerk_class.f90)' + + ! Get current particle state + state = p + + if (p % fate == LEAK_FATE) return + + ! 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 + scoreVal = self % response(1) % get(p, xsData) * p % w * flx + ! Deal with infinite cross sections - may not be the right solution for generality + if (scoreVal /= scoreVal) then + scoreVal = p % w + end if + call mem % score(scoreVal, adrr + MATERIAL_ENERGY) + + end subroutine reportHist + + !! + !! Tally the weight of a particle surviving the time step, corresponding to radiation energy in + !! that instant as a + !! + subroutine reportTrans(self, p, xsData, mem) + class(energyWeightClerk), 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 = 'reportTrans (energyWeightClerk_class.f90)' + + ! Get current particle state + state = p + + ! Consider only radiation particles (those surviving timestep) + if (p % fate /= AGED_FATE) return + if (p % type /= P_PHOTON) return + + ! 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 + scoreVal = self % response(1) % get(p, xsData) * p % w * flx + ! Deal with infinite cross sections - may not be the right solution for generality + if (scoreVal /= scoreVal) then + scoreVal = p % w + end if + call mem % score(scoreVal, adrr + RADIATION_ENERGY) + + end subroutine reportTrans + + !! + !! Display convergance progress on the console + !! + !! See tallyClerk_inter for details + !! + subroutine display(self, mem) + class(energyWeightClerk), intent(in) :: self + type(scoreMemory), intent(in) :: mem + + print *, 'energyWeightClerk 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(energyWeightClerk), 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 = [self % width, self % map % binArrayShape()] + else + resArrayShape = [self % width] + 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 + !! + !! See tallyClerk_inter for details + !! + pure subroutine getResult(self, res, mem) + class(energyWeightClerk), intent(in) :: self + class(tallyResult), allocatable, intent(inout) :: res + type(scoreMemory), intent(in) :: mem + real(defReal), dimension(:), allocatable :: mat, rad + integer(shortInt) :: i, N + + N = self % getSize() / 2 + allocate(mat(N)) + allocate(rad(N)) + + ! Get result value for each material + do i = 1, N + call mem % getResult(mat(i), self % getMemAddress()+2*i-2) + call mem % getResult(rad(i), self % getMemAddress()+2*i-1) + end do + + allocate(res, source = energyWeightClerkResult(mat,rad)) + + end subroutine getResult + +end module energyWeightClerk_class diff --git a/Tallies/TallyClerks/tallyClerkFactory_func.f90 b/Tallies/TallyClerks/tallyClerkFactory_func.f90 index f1277fa47..71a2003bc 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 energyWeightClerk_class, only : energyWeightClerk use mgXsClerk_class, only : mgXsClerk implicit none @@ -37,6 +38,7 @@ module tallyClerkFactory_func 'shannonEntropyClerk ',& 'centreOfMassClerk ',& 'dancoffBellClerk ',& + 'energyWeightClerk ',& 'mgXsClerk '] contains @@ -87,10 +89,13 @@ subroutine new_tallyClerk(new, dict, name) case('centreOfMassClerk') allocate(centreOfMassClerk :: new) + case('energyWeightClerk') + allocate(energyWeightClerk :: new) + case('mgXsClerk') allocate(mgXsClerk :: new) - case default + case default print *, AVALIBLE_tallyClerks call fatalError(Here, 'Unrecognised type of tallyClerk: ' // trim(type)) diff --git a/Tallies/TallyClerks/tallyClerk_inter.f90 b/Tallies/TallyClerks/tallyClerk_inter.f90 index 602fea2cb..b4ce4b642 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..87980261e 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 @@ -17,7 +19,7 @@ module weightResponse_class !! !! tallyResponse for scoring particle weights - !! Currently supports neutrons only + !! Currently supports neutrons and IMC only !! !! Interface: !! tallyResponse interface @@ -63,24 +65,22 @@ end subroutine init !! See tallyResponse_inter for details !! !! Errors: - !! Return ZERO if particle is not a Neutron + !! Return ZERO if particle is not a Neutron or IMC !! function get(self, p, xsData) result(val) class(weightResponse), intent(in) :: self 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 0e35adfdf..665856b21 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 @@ -749,6 +751,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..da7047c9d 100644 --- a/TransportOperator/CMakeLists.txt +++ b/TransportOperator/CMakeLists.txt @@ -4,4 +4,7 @@ add_sources(./transportOperator_inter.f90 ./transportOperatorDT_class.f90 # ./transportOperatorDynamicDT_class.f90 ./transportOperatorST_class.f90 - ./transportOperatorHT_class.f90) + ./transportOperatorHT_class.f90 + ./transportOperatorTime_class.f90 + ./transportOperatorGeomHT_class.f90 + ./virtualMat_class.f90) diff --git a/TransportOperator/transportOperatorFactory_func.f90 b/TransportOperator/transportOperatorFactory_func.f90 index a308ee01e..57989d293 100644 --- a/TransportOperator/transportOperatorFactory_func.f90 +++ b/TransportOperator/transportOperatorFactory_func.f90 @@ -12,7 +12,8 @@ module transportOperatorFactory_func use transportOperatorST_class, only : transportOperatorST use transportOperatorDT_class, only : transportOperatorDT use transportOperatorHT_class, only : transportOperatorHT - !use transportOperatorDynamicDT_class, only : transportOperatorDynamicDT + use transportOperatorTime_class, only : transportOperatorTime + use transportOperatorGeomHT_class, only : transportOperatorGeomHT implicit none private @@ -21,9 +22,11 @@ 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 ', & + 'transportOperatorTime ', & + 'transportOperatorGeomHT'] public :: new_transportOperator @@ -44,7 +47,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 select case(type) case('transportOperatorST') allocate( transportOperatorST :: new) @@ -55,6 +58,12 @@ subroutine new_transportOperator(new, dict) case('transportOperatorHT') allocate( transportOperatorHT :: new) + case('transportOperatorTime') + allocate( transportOperatorTime :: new) + + case('transportOperatorGeomHT') + allocate( transportOperatorGeomHT :: new) + case default print *, AVALIBLE_transportOps call fatalError(Here, 'Unrecognised type of transportOperator: ' // trim(type)) diff --git a/TransportOperator/transportOperatorGeomHT_class.f90 b/TransportOperator/transportOperatorGeomHT_class.f90 new file mode 100644 index 000000000..63fe2ad4a --- /dev/null +++ b/TransportOperator/transportOperatorGeomHT_class.f90 @@ -0,0 +1,377 @@ +!! +!! Hybrid transport operator using two geometries to switch between delta tracking and surface tracking +!! +module transportOperatorGeomHT_class + use numPrecision + use universalVariables + + use genericProcedures, only : fatalError, numToChar + use particle_class, only : particle, P_PHOTON, P_MATERIAL + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary + use RNG_class, only : RNG + + ! Superclass + use transportOperator_inter, only : transportOperator, init_super => init + + ! Geometry interfaces + use geometry_inter, only : geometry + use geometryGrid_class, only : geometryGrid + use geometryReg_mod, only : gr_geomPtr => geomPtr, & + gr_geomIdx => geomIdx + use geometryFactory_func, only : new_geometry + use coord_class, only : coordList + + ! Tally interface + use tallyCodes + use tallyAdmin_class, only : tallyAdmin + + ! Nuclear data interfaces + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase, mgIMCDatabase_CptrCast + use nuclearDataReg_mod, only : ndReg_get => get + + use virtualMat_class, only : virtualMat + use simulationTime_class + + implicit none + private + + !! + !! Transport operator that moves a particle using hybrid tracking, up to a time boundary. + !! + !! Overlays a second geometry (currently only able to be geometryGrid_class) onto the primary + !! geometry. This second geometry is filled with fake materials (see virtualMat_class) to give + !! information about the primary geometry. A particle begins with delta tracking, using the + !! majorant of only the current cell in the upper (fake) geometry. If it would cross into a new + !! upper cell, it is moved in the upper geometry and the majorant is used accordingly. If the + !! acceptance probability is ever below the user-supplied cutoff (actually 1 - cutoff), standard + !! surface tracking is used in the lower geometry, for the rest of the particle's motion. + !! + !! In order to determine which real mats are in which virtual cell, input searchN = (Nx, Ny, Nz) + !! is provided, and Nx*Ny*Nz points are generated and used to map between the two geometries. If + !! searchN is too small then materials may be missed! (fatalError will be called if local opacity + !! is greater than majorant opacity for any particle, which will be the result of missed mats). + !! If searchN = (N) is given (size-1 array) then Nx = Ny = Nz = N. + !! + !! Sample Dictionary Input: + !! transportOperator { + !! type transportOperatorGeomHT; + !! cutoff 0.5; + !! geometry { type geometryGrid; + !! bounds (-2.4 -0.5 -0.5 2.4 0.5 0.5); => Should match primary geometry + !! boundary (1 1 1 1 1 1); => Should match primary geometry + !! gridOnly y; => Stops materials being overwritten + !! dimensions (25 1 1); } => Number of tracking cells + !! searchN (100); + !! } + !! + !! TODO There is a bug somewhere, potentially in this module, causing segmentation faults in + !! certain situations, but I haven't been able to figure out where and why. + !! See InputFiles/IMC/DensmoreMF/densmoreMid for the input that leads to this fault. + !! + type, public, extends(transportOperator) :: transportOperatorGeomHT + real(defReal) :: deltaT + real(defReal) :: cutoff + integer(shortInt) :: method + class(geometry), pointer :: upperGeom + integer(shortInt) :: upperGeomIdx + integer(shortInt) :: thisTimeStep + class(virtualMat), dimension(:), pointer :: virtualMats + contains + procedure :: transit => multiGeomTracking + procedure :: init + procedure, private :: surfaceTracking + procedure, private :: deltaTracking + end type transportOperatorGeomHT + +contains + + !! + !! Update virtual materials if needed, and then begin particle journey with delta tracking. + !! + subroutine multiGeomTracking(self, p, tally, thisCycle, nextCycle) + class(transportOperatorGeomHT), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + class(particleDungeon), intent(inout) :: thisCycle + class(particleDungeon), intent(inout) :: nextCycle + integer(shortInt) :: i + character(100), parameter :: Here = 'timeTracking (transportOperatorGeomHT_class.f90)' + + ! Material transform not included in this class to avoid complicating any further. Can be + ! essentially copied and pasted from transportOperatorTime_class if desired to use ISMC here + if (p % type == P_MATERIAL) call fatalError(Here, 'No support for ISMC in this transOp') + + ! Update majorants if required - this would be better done at the end of time step in PP to + ! avoid check for each particle and repetion in parallel but I wanted to keep this class + ! self-contained for now + if (self % thisTimeStep /= thisStep()) then + self % thisTimeStep = thisStep() + do i = 1, size(self % virtualMats) + call self % virtualMats (i) % updateMajorant() + end do + end if + + call self % deltaTracking(p) + + ! 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 multiGeomTracking + + !! + !! Perform delta tracking + !! + !! Note that the method used of dealing with upper geometry cell changes does lead to potential + !! inaccuracy if a particle can change upperGeom cells and then return to the original one in the + !! same move, for example if upperGeom is curved or if there is a reflection. For now this class + !! requires upperGeom to be geometryGrid_class (no curves) and this sort of reflection should + !! never happen for any vaguely sensible geometry/material choices. One could instead explicitly + !! calculate distance to a new upperGeom cell and compare to dColl and dTime, at the cost of + !! efficiency. + !! + subroutine deltaTracking(self, p) + class(transportOperatorGeomHT), intent(inout) :: self + class(particle), intent(inout) :: p + class(coordList), allocatable :: coords + real(defReal) :: dTime, dColl, sigmaT, majorant_inv, dist, ratio + integer(shortInt) :: virtualMatIdx, testMat, uniqueID, event + character(100), parameter :: Here = 'deltaTracking (transportOperatorGeomHT_class.f90)' + + ! Get index of virtual material and majorant inverse + call self % upperGeom % whatIsAt(virtualMatIdx, uniqueID, p % coords % lvl(1) % r) + majorant_inv = self % virtualMats(virtualMatIdx) % majorant_inv(p % G) + + ! Get initial local opacity + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + + DTLoop:do + + ! Switch to surface tracking if delta tracking is unsuitable + ratio = sigmaT * majorant_inv + if (ratio > ONE) call fatalError(Here, 'Local opacity greater than majorant') + if (ratio <= self % cutoff .or. majorant_inv == ZERO) then + call self % surfaceTracking(p) + return + end if + + ! Find distance to time boundary + dTime = lightSpeed * (p % timeMax - p % time) + + ! Sample distance to collision + dColl = -log(p % pRNG % get()) * majorant_inv + + ! Move copy of coords by minimum distance without considering surface crossings + coords = p % coords + dist = min(dColl, dTime) + call self % geom % teleport(coords, dist) + + ! Check for particle leakage + if (coords % matIdx == OUTSIDE_FILL) then + p % coords = coords + return + end if + + ! Check for change of upper geometry + call self % upperGeom % whatIsAt(testMat, uniqueID, coords % lvl(1) % r) + if (testMat /= virtualMatIdx) then + ! Move would take particle to a new cell + call self % upperGeom % move(p % coords, dist, event) + ! Get new majorant (particle already placed in upper geometry) + virtualMatIdx = p % matIdx() + majorant_inv = self % virtualMats(virtualMatIdx) % majorant_inv(p % G) + ! Update particle time and place back in lower geometry + p % time = p % time + dist / lightSpeed + call self % geom % placeCoord(p % coords) + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + ! Return to start of loop + cycle DTLoop + end if + + ! Move accepted, move p + p % coords = coords + + ! Update particle time + p % time = p % time + dist / lightSpeed + + ! Act based on distance moved + if (dist == dTime) then + ! Update particle fate and exit + p % fate = AGED_FATE + p % time = p % timeMax + exit DTLoop + + else if (dist == dColl) then + ! Get local opacity and check for real or virtual collision + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + if (p % pRNG % get() < sigmaT * majorant_inv) exit DTLoop + + else + call fatalError(Here, 'Peculiar distance moved') + end if + + end do DTLoop + + end subroutine deltaTracking + + !! + !! Perform standard surface tracking using only the lower geometry. + !! Once in this loop, delta tracking is not used again. + !! + subroutine surfaceTracking(self, p) + class(transportOperatorGeomHT), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime, dColl, dist, sigmaT + integer(shortInt) :: event + character(100), parameter :: Here = 'surfaceTracking (transportOperatorGeomHT_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 + + ! Ensure particle does not remain exactly on a boundary if dColl is close to 0 + if (event == CROSS_EV .and. dColl < SURF_TOL) then + dColl = SURF_TOL + end if + + ! 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 + p % fate = AGED_FATE + p % time = p % timeMax + exit STLoop + + else if (dist == dColl) then + ! Collision + exit STLoop + + end if + + if (event == COLL_EV) call fatalError(Here, 'Move outcome should be CROSS_EV or BOUNDARY_EV') + + end do STLoop + + if (event /= COLL_EV) call fatalError(Here, 'Move outcome should be COLL_EV') + + end subroutine surfaceTracking + + !! + !! 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(transportOperatorGeomHT), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(nameLen) :: geomName + class(dictionary),pointer :: tempdict + class(geometry), pointer :: upperGeom + integer(shortInt), dimension(:), allocatable :: dimensions, searchN + integer(shortInt), dimension(3) :: searchN3 + integer(shortInt) :: N, i, j, k + integer(shortInt) :: realMatIdx, virtualMatIdx, uniqueID + real(defReal), dimension(6) :: bounds + real(defReal), dimension(3) :: corner, extent, searchRes, r + character(100), parameter :: Here = "init (transportOperatorGeomHT_class.f90)" + + ! Initialise superclass + call init_super(self, dict) + + ! Get cutoff value + call dict % getOrDefault(self % cutoff, 'cutoff', 0.3_defReal) + ! Flip to be consistent with transportOperatorHT_class + self % cutoff = ONE - self % cutoff + + ! Build upper level geometry + geomName = 'surfaceGeom' + tempDict => dict % getDictPtr('geometry') + call new_geometry(tempDict, geomName) + self % upperGeomIdx = gr_geomIdx(geomName) + upperGeom => gr_geomPtr(self % upperGeomIdx) + self % upperGeom => upperGeom + + ! Provide access to lower (standard) geometry + ! TODO: Really geometry and xsData should be inputs to init, but would need to change all + ! other transport operators + self % geom => gr_geomPtr(1) + self % xsData => ndReg_get(P_PHOTON_MG) + + ! For now limited to grid geometry + select type(upperGeom) + class is(geometryGrid) + ! Get some basic geometry info + corner = upperGeom % corner + class default + call fatalError(Here, 'Geometry class should be geometryGrid') + end select + + ! Initialise a virtual material object for each cell + call tempDict % get(dimensions,'dimensions') + N = dimensions(1) * dimensions(2) * dimensions(3) + allocate(self % virtualMats(N)) + do i = 1, N + call self % virtualMats(i) % init(self % xsData) + end do + + ! Get resolution to search through grid + call dict % get(searchN, 'searchN') + if (size(searchN) == 3) then + searchN3 = searchN + else if (size(searchN) == 1) then + do i = 1, 3 + searchN3(i) = searchN(1) + end do + else + call fatalError(Here, 'searchN must be of size 1 or 3') + end if + if (any(searchN3 < 1)) call fatalError(Here, 'Invalid searchN') + + ! Search grid to assign real materials to virtual materials + bounds = upperGeom % bounds() + do i = 1, 3 + extent(i) = bounds(i+3) - bounds(i) + end do + searchRes = extent / (searchN3 + 1) + do i = 1, searchN3(1) + do j = 1, searchN3(2) + do k = 1, searchN3(3) + ! Find matIdx at search location + r = corner + [i, j, k] * searchRes + call self % geom % whatIsAt(realMatIdx, uniqueID, r) + call upperGeom % whatIsAt(virtualMatIdx, uniqueID, r) + call self % virtualMats(virtualMatIdx) % addRealMat(realMatIdx) + end do + end do + end do + + do i = 1, size(self % virtualMats) + call self % virtualMats (i) % updateMajorant() + end do + self % thisTimeStep = thisStep() + + end subroutine init + +end module transportOperatorGeomHT_class diff --git a/TransportOperator/transportOperatorIMC_class.f90 b/TransportOperator/transportOperatorIMC_class.f90 new file mode 100644 index 000000000..95ce041b8 --- /dev/null +++ b/TransportOperator/transportOperatorIMC_class.f90 @@ -0,0 +1,495 @@ +!! +!! Transport operator for implicit Monte Carlo scheme using delta tracking +!! +module transportOperatorIMC_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 + + ! Geometry interfaces + use geometry_inter, only : geometry + + ! Tally interface + use tallyCodes + use tallyAdmin_class, only : tallyAdmin + + ! Nuclear data interfaces + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCMaterial_inter, only : mgIMCMaterial, mgIMCMaterial_CptrCast + use materialMenu_mod, only : mm_nMat => nMat + + implicit none + private + + !! + !! Transport operator that moves a particle with IMC tracking + !! + type, public, extends(transportOperator) :: transportOperatorIMC + class(mgIMCMaterial), pointer, public :: mat => null() + real(defReal) :: majorant_inv + real(defReal) :: deltaT + real(defReal) :: cutoff + real(defReal), dimension(:), allocatable :: matMajs + real(defReal), dimension(:), allocatable :: sigmaLocal + integer(shortInt), dimension(:,:), allocatable :: matConnections + integer(shortInt) :: majMapN = 0 + real(defReal), dimension(3) :: top = ZERO + real(defReal), dimension(3) :: bottom = ZERO + integer(shortInt) :: pSteps + contains + procedure :: transit => imcTracking + procedure :: init + procedure :: buildMajMap + procedure :: updateMajorants + procedure, private :: materialTransform + procedure, private :: surfaceTracking + procedure, private :: deltaTracking + procedure, private :: simpleParticle + end type transportOperatorIMC + +contains + + subroutine imcTracking(self, p, tally, thisCycle, nextCycle) + class(transportOperatorIMC), 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 = 'IMCTracking (transportOperatorIMC_class.f90)' + + ! Deal with material particles, only relevant for ISMC + if(p % getType() == P_MATERIAL_MG) then + call self % materialTransform(p, tally) + if(p % fate == AGED_FATE) return + end if + + ! Get majorant for particle + if (allocated(self % matMajs)) then + self % majorant_inv = ONE / self % matMajs(p % matIdx()) + else + self % majorant_inv = ONE / self % xsData % getMajorantXS(p) + end if + + ! Check for errors + if (p % getType() /= P_PHOTON_MG) call fatalError(Here, 'Particle is not MG Photon') + if (p % time /= p % time) call fatalError(Here, 'Particle time is NaN') + + ! Obtain sigmaT + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + + if (sigmaT*self % majorant_inv > 1) call fatalError(Here, 'Sigma greater than majorant.& + & MajorantMap settings may have been chosen poorly.') + + ! Decide whether to use delta tracking or surface tracking + ! Vastly different opacities make delta tracking infeasable + if(sigmaT * self % majorant_inv > ONE - self % cutoff .or. self % cutoff == ONE) then + ! Delta tracking + call self % deltaTracking(p) + else + ! Surface tracking + call self % surfaceTracking(p) + end if + + ! Check for particle leakage + if (p % matIdx() == OUTSIDE_FILL) then + p % fate = LEAK_FATE + p % isDead = .true. + return + end if + + call tally % reportTrans(p) + + end subroutine imcTracking + + !! + !! Perform surface tracking + !! + subroutine surfaceTracking(self, p) + class(transportOperatorIMC), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime + real(defReal) :: dColl + real(defReal) :: dist, sigmaT + integer(shortInt) :: event + character(100), parameter :: Here = 'surfaceTracking (transportOperatorIMC_class.f90)' + + STLoop:do + + ! Find distance to time boundary + dTime = lightSpeed * (p % timeMax - p % time) + + ! Sample distance to collision + if (p % matIdx() == VOID_MAT) then + dColl = INF + else + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + dColl = -log( p % pRNG % get() ) / sigmaT + end if + + ! 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 + if (abs(p % time - p % timeMax)>0.000001) call fatalError(Here, 'Particle time is somehow incorrect') + 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 + !! + subroutine deltaTracking(self, p) + class(transportOperatorIMC), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime + real(defReal) :: dColl + real(defReal) :: sigmaT + character(100), parameter :: Here = 'deltaTracking (transportOperatorIMC_class.f90)' + + DTLoop:do + + ! Find distance to time boundary + dTime = lightSpeed * (p % timeMax - p % time) + + ! Sample distance to collision + if (p % matIdx() == VOID_MAT) then + dColl = INF + else + dColl = -log( p % pRNG % get() ) * self % majorant_inv + end if + + ! If dTime < dColl, move to end of time step location + if (dTime < dColl) then + call self % geom % teleport(p % coords, dTime) + p % fate = AGED_FATE + p % time = p % timeMax + exit DTLoop + end if + + ! Otherwise, move to potential collision location + call self % geom % teleport(p % coords, dColl) + p % time = p % time + dColl / lightSpeed + + ! Check for particle leakage + if (p % matIdx() == OUTSIDE_FILL) return + + ! Obtain local cross-section + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) + + ! Roll RNG to determine if the collision is real or virtual + ! Exit the loop if the collision is real + if (p % pRNG % get() < sigmaT * self % majorant_inv) exit DTLoop + + ! Protect against infinite loop + if (sigmaT * self % majorant_inv == 0) call fatalError(Here, '100 % virtual collision chance,& + & potentially infinite loop') + + ! Switch to ST if particle moves into a region where delta tracking is no longer feasible + if (sigmaT * self % majorant_inv < ONE - self % cutoff) then + call self % surfaceTracking(p) + return + end if + + end do DTLoop + + end subroutine deltaTracking + + !! + !! Transform material particles into radiation photons with + !! probability per unit time of c*sigma_a*fleck*eta + !! + !! Used only for ISMC, not for standard IMC + !! + subroutine materialTransform(self, p, tally) + class(transportOperatorIMC), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + real(defReal) :: sigmaT, fleck, eta, mu, phi + real(defReal), dimension(3) :: dir + character(100), parameter :: Here = 'materialTransform (transportOperatorIMC_class.f90)' + + ! 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") + + sigmaT = self % xsData % getTransMatXS(p, p % matIdx()) !! Should be sigma_a, may need changing when sorting out cross-sections + fleck = self % mat % getFleck() + eta = self % mat % getEta() + + ! Sample time to transform into radiation photon + p % time = p % time - log(p % pRNG % get()) / (sigmaT*fleck*eta*lightSpeed) + + ! Deal with eta = 0 causing NaN + if (p % time /= p % time) p % time = p % time -log(p % pRNG % get()) / (1.400*fleck*lightSpeed) + + ! Exit loop if particle remains material until end of time step + if (p % time >= p % timeMax) then + p % fate = AGED_FATE + p % time = p % timeMax + ! Tally energy for next temperature calculation + call tally % reportHist(p) + else + p % type = P_PHOTON + ! Resample direction + 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) + call p % point(dir) + end if + + end subroutine materialTransform + + !! + !! Generate particles and follow them to attempt to reduce majorant opacity seen by each cell. + !! + !! Particles sampled uniformly within geometry, with random directions, then incrementally moved + !! up to distance dTime. Moving increments are determined by the number of steps to be taken, + !! either given directly in input or calculated from dTime and a given cell lengthscale. + !! + !! Majorant of each cell is set to be the maximum opacity seen by any particle originating in + !! that cell. + !! + subroutine buildMajMap(self, rand, xsData) + class(transportOperatorIMC), intent(inout) :: self + class(RNG), intent(inout) :: rand + class(nuclearDatabase), intent(in), pointer :: xsData + type(particle) :: p + integer(shortInt) :: i, j, matIdx + real(defReal) :: dist + + ! Check that subroutine should be called + if (.not. allocated(self % matMajs)) return + + ! Point to nuclear data, as otherwise cannot access until after first particle transport + self % xsData => xsData + + ! Reset array + self % matConnections = 0 + + ! Calculate distance increments + dist = self % deltaT * lightSpeed / self % pSteps + + do i = 1, self % majMapN + + ! Sample particle + call self % simpleParticle(p, rand) + matIdx = p % matIdx() + + ! Incrementally transport particle up to a distance dTime + do j = 1, self % pSteps + + call self % geom % teleport(p % coords, dist) + if (p % matIdx() == VOID_MAT .or. p % matIdx() == OUTSIDE_MAT) exit + + ! Update matConnections to signify a connection between starting mat and new mat + self % matConnections(matIdx, p % matIdx()) = 1 + self % matConnections(p % matIdx(), matIdx) = 1 + + end do + + end do + + end subroutine buildMajMap + + !! + !! Update majorants for each region using material connections build up in buildMajMap subroutine + !! + subroutine updateMajorants(self, rand) + class(transportOperatorIMC), intent(inout) :: self + class(RNG), intent(inout) :: rand + integer(shortInt) :: i, G, nMats + character(100), parameter :: Here = 'updateMajorants (transportOperatorIMC_class.f90)' + + ! Check that subroutine should be called + if (.not. allocated(self % matMajs)) return + + nMats = mm_nMat() + G = 1 ! Can easily be extended to multiple groups later + + ! First, update array of local opacities + do i = 1, nMats + ! Get and verify material pointer for material i + self % mat => mgIMCMaterial_CptrCast(self % xsData % getMaterial(i)) + if(.not.associated(self % mat)) call fatalError(Here, "Failed to get MG IMC Material") + + ! Store opacity + self % sigmaLocal(i) = self % mat % getTotalXS(G, rand) + end do + + ! Now update majorants for each material + do i = 1, nMats + self % matMajs(i) = maxval(self % sigmaLocal * self % matConnections(i, 1:nMats)) + end do + + !print *, 'Local opacities:' + !print *, self % sigmaLocal + + !print *, 'New majorants:' + !print *, self % matMajs + + end subroutine updateMajorants + + !! + !! Sample position for buildMajMap subroutine (see above) + !! Attach only necessary properties to particle: + !! + !! - Position, sampled uniformly within geometry + !! - Direction, uniformly from unit sphere + !! - Type = P_PHOTON + !! - Group = 1, can easily extend to work with multiple groups another time + !! + subroutine simpleParticle(self, p, rand) + class(transportOperatorIMC), intent(inout) :: self + type(particle), intent(inout) :: p + class(RNG), intent(inout) :: rand + real(defReal) :: mu, phi + real(defReal), dimension(3) :: r, dir, rand3 + integer(shortInt) :: matIdx, uniqueID, loops + character(100), parameter :: Here = 'simpleParticle (transportOperatorIMC.f90)' + + ! Sample points randomly within geometry until valid material is found + loops = 0 + positionSample:do + ! Protect against infinite loop + loops = loops + 1 + if (loops >= 500) call fatalError(Here, '500 particles sampled in void or outside geometry') + + ! Sample position + rand3(1) = rand % get() + rand3(2) = rand % get() + rand3(3) = rand % get() + r = (self % top - self % bottom) * rand3 + self % bottom + + ! Find material under position + call self % geom % whatIsAt(matIdx, uniqueID, r) + + ! Reject if there is no material + if (matIdx == VOID_MAT .or. matIdx == OUTSIDE_MAT) cycle positionSample + + call p % coords % assignPosition(r) + exit positionSample + + end do positionSample + + ! 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) + call p % coords % assignDirection(dir) + + p % type = P_PHOTON + p % G = 1 + p % isMG = .true. + + call self % geom % placeCoord(p % coords) + + end subroutine simpleParticle + + !! + !! Get transport settings + !! + !! Sample dictionary input: + !! + !! transportOperator { + !! type transportOperatorIMC; + !! cutoff 0.5; + !! majMap { + !! nParticles 500; + !! pSteps 10; + !! } + !! } + !! + !! Cutoff of 1 gives exclusively delta tracking, cutoff of 0 gives exclusively surface tracking + !! + !! As an alternative to 'pSteps' can specify 'lengthScale' and then steps is calculated + !! automatically as pSteps = c*dt/lengthScale + !! + subroutine init(self, dict, geom) + class(transportOperatorIMC), intent(inout) :: self + class(dictionary), intent(in) :: dict + class(geometry), pointer, intent(in), optional :: geom + class(dictionary), pointer :: tempDict + integer(shortInt) :: nMats + real(defReal), dimension(6) :: bounds + real(defReal) :: lengthScale + character(100), parameter :: Here = 'init (transportOperatorIMC.f90)' + + ! Initialise superclass + call init_super(self, dict) + + self % geom => geom + + ! Get timestep size + call dict % get(self % deltaT, 'deltaT') + + ! Get cutoff value + call dict % getOrDefault(self % cutoff, 'cutoff', 0.7_defReal) + + ! Preparation for majorant reduction subroutine + if (dict % isPresent('majMap')) then + + ! Get settings + tempDict => dict % getDictPtr('majMap') + call tempDict % get(self % majMapN, 'nParticles') + + if (tempDict % isPresent('pSteps')) then + call tempDict % get(self % pSteps, 'pSteps') + else + call tempDict % get(lengthScale, 'lengthScale') + self % pSteps = ceiling(lightSpeed*self % deltaT/lengthScale) + end if + + nMats = mm_nMat() + + ! Allocate arrays + allocate(self % matMajs(nMats)) + allocate(self % sigmaLocal(nMats)) + allocate(self % matConnections(nMats, nMats)) + self % matMajs = 0 + self % sigmaLocal = 0 + + ! Set bounding region for particle sourcing + bounds = self % geom % bounds() + self % bottom = bounds(1:3) + self % top = bounds(4:6) + end if + + end subroutine init + +end module transportOperatorIMC_class diff --git a/TransportOperator/transportOperatorTime_class.f90 b/TransportOperator/transportOperatorTime_class.f90 new file mode 100644 index 000000000..cda201b59 --- /dev/null +++ b/TransportOperator/transportOperatorTime_class.f90 @@ -0,0 +1,185 @@ +!! +!! Surface tracking transport operator for time-dependent problems +!! +module transportOperatorTime_class + use numPrecision + use universalVariables + + use genericProcedures, only : fatalError, numToChar + use particle_class, only : particle, P_PHOTON, P_MATERIAL + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary + use rng_class, only : rng + + ! Superclass + use transportOperator_inter, only : transportOperator + + ! Geometry interfaces + use geometry_inter, only : geometry + + ! Tally interface + use tallyCodes + use tallyAdmin_class, only : tallyAdmin + + ! Nuclear data interfaces + use nuclearDatabase_inter, only : nuclearDatabase + use mgIMCDatabase_inter, only : mgIMCDatabase, mgIMCDatabase_CptrCast + + implicit none + private + + !! + !! Transport operator that moves a particle using surface tracking, up to a time boundary + !! + type, public, extends(transportOperator) :: transportOperatorTime + contains + procedure :: transit => timeTracking + procedure, private :: surfaceTracking + procedure, private :: materialTransform + end type transportOperatorTime + +contains + + subroutine timeTracking(self, p, tally, thisCycle, nextCycle) + class(transportOperatorTime), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + class(particleDungeon), intent(inout) :: thisCycle + class(particleDungeon), intent(inout) :: nextCycle + character(100), parameter :: Here = 'timeTracking (transportOperatorTime_class.f90)' + + ! Transform material particles into photons - for use in ISMC + if (p % type == P_MATERIAL) then + call self % materialTransform(p, tally) + ! Exit at time boundary + if (p % fate == AGED_FATE) return + end if + + ! Perform tracking + call self % surfaceTracking(p) + + ! 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(transportOperatorTime), intent(inout) :: self + class(particle), intent(inout) :: p + real(defReal) :: dTime, dColl, dist, sigmaT + integer(shortInt) :: event + character(100), parameter :: Here = 'surfaceTracking (transportOperatorTime_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 + + ! Ensure particle does not remain exactly on a boundary if dColl is close to 0 + if (event == CROSS_EV .and. dColl < SURF_TOL) then + dColl = SURF_TOL + end if + + ! 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 + p % fate = AGED_FATE + p % time = p % timeMax + exit STLoop + + else if (dist == dColl) then + ! Collision + exit STLoop + + end if + + if (event == COLL_EV) call fatalError(Here, 'Move outcome should be CROSS_EV or BOUNDARY_EV') + + end do STLoop + + if (event /= COLL_EV) call fatalError(Here, 'Move outcome should be COLL_EV') + + end subroutine surfaceTracking + + !! + !! Determine when a material particle will transform into a photon for ISMC calculations + !! + !! Args: + !! p [inout] -> material particle to be transformed + !! tally [inout] -> tally to keep track of material particles surviving time step + !! + subroutine materialTransform(self, p, tally) + class(transportOperatorTime), intent(inout) :: self + class(particle), intent(inout) :: p + type(tallyAdmin), intent(inout) :: tally + real(defReal) :: transformTime, mu, phi + real(defReal), dimension(3) :: dir + class(mgIMCDatabase), pointer :: nucData + integer(shortInt) :: matIdx, uniqueID + character(100), parameter :: Here = 'materialTransform (transportOperatorIMC_class.f90)' + + ! Get pointer to nuclear database + nucData => mgIMCDatabase_CptrCast(self % xsData) + if (.not. associated(nucData)) call fatalError(Here, 'Unable to find mgIMCDatabase') + + ! Material particles can occasionally have coords placed in void if within SURF_TOL of boundary + matIdx = p % matIdx() + ! If so, get matIdx based on exact position (no adjustment for surface tol) + if (matIdx == VOID_MAT .or. matIdx == OUTSIDE_MAT) then + call self % geom % whatIsAt(matIdx, uniqueID, p % coords % lvl(1) % r, [ZERO,ZERO,ZERO]) + end if + ! If still in invalid region, call fatalError + if (matIdx == 0) call fatalError(Here, 'Outside material particle') + if (matIdx == VOID_MAT) call fatalError(Here, 'Void material particle') + + ! Sample time until emission + transformTime = nucData % sampleTransformTime(matIdx, p % pRNG) + p % time = min(p % timeMax, p % time + transformTime) + + ! Exit loop if particle remains material until end of time step + if (p % time == p % timeMax) then + p % fate = AGED_FATE + ! Tally energy for next temperature calculation + call tally % reportHist(p) + + ! Transform into photon + else + p % type = P_PHOTON + ! Resample direction + 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) + call p % point(dir) + ! Resample energy + p % G = nucData % sampleEnergyGroup(matIdx, p % pRNG) + end if + + end subroutine materialTransform + +end module transportOperatorTime_class diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index 8d0e30e68..aaeba6226 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -7,6 +7,7 @@ module transportOperator_inter use particle_class, only : particle use particleDungeon_class, only : particleDungeon use dictionary_class, only : dictionary + use RNG_class, only : RNG ! Geometry interfaces use geometryReg_mod, only : gr_geomPtr => geomPtr @@ -19,8 +20,6 @@ module transportOperator_inter use nuclearDataReg_mod, only : ndReg_get => get use nuclearDatabase_inter, only : nuclearDatabase - - implicit none private @@ -121,8 +120,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 @@ -139,5 +138,4 @@ elemental subroutine kill(self) end subroutine kill - end module transportOperator_inter diff --git a/TransportOperator/virtualMat_class.f90 b/TransportOperator/virtualMat_class.f90 new file mode 100644 index 000000000..06111eb1d --- /dev/null +++ b/TransportOperator/virtualMat_class.f90 @@ -0,0 +1,111 @@ +module virtualMat_class + + use numPrecision + use universalVariables + + use genericProcedures, only : fatalError, numToChar + + use dynArray_class, only : dynIntArray + + use nuclearDatabase_inter, only : nuclearDatabase + + use particle_class, only : particle, P_PHOTON, P_MATERIAL + + use materialEquations, only : mgEnergyGrid + + implicit None + private + + !! + !! Fake material class used for multiGeom delta tracking. Fills upper geometry cells, and stores + !! information about the real materials overlapped in the lower geometry, specifically which + !! materials these are and the majorant_inv within the upper geometry cell. + !! + !! See transportOperatorGeomHT_class.f90 for example of use. + !! + type, public :: virtualMat + type(dynIntArray) :: realMats + real(defReal), dimension(:), allocatable :: majorant_inv + class(nuclearDatabase), pointer :: xsData + contains + procedure :: init + procedure :: addRealMat + procedure :: updateMajorant + end type virtualMat + +contains + + !! + !! Allocate space for majorants + !! + subroutine init(self, nucData) + class(virtualMat), intent(inout) :: self + class(nuclearDatabase), pointer, intent(in) :: nucData + integer(shortInt) :: nG + + self % xsData => nucData + + ! Allocate space for majorants + ! TODO: Currently this method of obtaining nG is specific to MG IMC + if(associated(mgEnergyGrid)) then + nG = mgEnergyGrid % getSize() + else + nG = 1 + end if + + allocate(self % majorant_inv(nG)) + + end subroutine + + !! + !! Add matIdx to the list of real materials in the underlying geometry that overlap with the + !! region occupied by this virtual material. + !! + subroutine addRealMat(self, matIdx) + class(virtualMat), intent(inout) :: self + integer(shortInt), intent(in) :: matIdx + + ! Add matIdx to array if not already present + if (.not. self % realMats % isPresent(matIdx)) then + call self % realMats % add(matIdx) + end if + + end subroutine addRealMat + + !! + !! Update the majorant_inv for each group + !! + subroutine updateMajorant(self) + class(virtualMat), intent(inout) :: self + integer(shortInt) :: i, G + real(defReal) :: majorant, sigma + class(particle), allocatable :: p + character(100), parameter :: Here = '' + + ! Create particle for call to obtain XS - no data needed except p % G + allocate(p) + p % isMG = .true. + + majorant = ZERO + + ! Find majorant opacity of virtual material for each group + do G = 1, size(self % majorant_inv) + p % G = G + do i = 1, self % realMats % getSize() + sigma = self % xsData % getTransMatXS(p, self % realMats % get(i)) + if (sigma > majorant) majorant = sigma + end do + + ! Avoid infinite maj_inv for virtualMat cells containing only void + if (majorant == ZERO) then + ! maj_inv = 0 forces particle to default to ST in transport operator to avoid issues + self % majorant_inv(G) = ZERO + else + self % majorant_inv(G) = ONE/majorant + end if + + end do + + end subroutine updateMajorant + +end module virtualMat_class