diff --git a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 index c954729b5..5c8c29208 100644 --- a/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 +++ b/CollisionOperator/CollisionProcessors/collisionProcessor_inter.f90 @@ -16,7 +16,7 @@ module collisionProcessor_inter private !! - !! Data package with all relevant data about the collision to move beetween customisable + !! Data package with all relevant data about the collision to move between customisable !! procedures !! type, public :: collisionData diff --git a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 index fd8ce8058..3a9d815d4 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEimp_class.f90 @@ -2,7 +2,7 @@ module neutronCEimp_class use numPrecision use endfConstants - use universalVariables, only : nameUFS, nameWW, REJECTED + use universalVariables, only : nameUFS, nameWW, REJECTED, kBoltzmannMev use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -241,7 +241,7 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle type(neutronMicroXSs) :: microXSs - real(defReal) :: r, denom, alphaXS, probAlpha + real(defReal) :: r, kT, denom, alphaXS, probAlpha character(100),parameter :: Here = 'sampleCollision (neutronCEimp_class.f90)' ! Verify that particle is CE neutron @@ -273,7 +273,7 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) if(.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial') ! Select collision nuclide - call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E) + call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E, p % T, p % rho) ! If nuclide was rejected in TMS loop return to tracking if (collDat % nucIdx == REJECTED) then @@ -285,7 +285,12 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) if (.not.associated(self % mat)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') ! Select Main reaction channel - call self % nuc % getMicroXSs(microXss, collDat % E, self % mat % kT, p % pRNG) + if (p % T <= ZERO) then + kT = self % mat % kT + else + kT = p % T * kBoltzmannMeV + end if + call self % nuc % getMicroXSs(microXss, collDat % E, kT, p % pRNG) r = p % pRNG % get() collDat % MT = microXss % invert(r) @@ -308,7 +313,7 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) integer(shortInt) :: n, i real(defReal) :: wgt, rand1, E_out, mu, phi, lambda real(defReal) :: sig_nufiss, sig_tot, k_eff, & - sig_scatter, totalElastic + sig_scatter, totalElastic, kT logical(defBool) :: fiss_and_implicit, keepDel character(100),parameter :: Here = 'implicit (neutronCEimp_class.f90)' @@ -323,7 +328,12 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) rand1 = p % pRNG % get() ! Random number to sample sites ! Retrieve cross section at the energy used for reaction sampling - call self % nuc % getMicroXSs(microXSs, collDat % E, self % mat % kT, p % pRNG) + if (p % T <= ZERO) then + kT = self % mat % kT + else + kT = p % T * kBoltzmannMeV + end if + call self % nuc % getMicroXSs(microXSs, collDat % E, kT, p % pRNG) sig_nufiss = microXSs % nuFission sig_tot = microXSs % total @@ -445,7 +455,7 @@ subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) real(defReal),dimension(3) :: r, dir, val integer(shortInt) :: n, i real(defReal) :: wgt, rand1, E_out, mu, phi, lambda - real(defReal) :: sig_nufiss, sig_fiss, k_eff, wD + real(defReal) :: sig_nufiss, sig_fiss, k_eff, kT, wD character(100),parameter :: Here = 'fission (neutronCEimp_class.f90)' if (.not.self % implicitSites) then @@ -456,7 +466,12 @@ subroutine fission(self, p, tally, collDat, thisCycle, nextCycle) rand1 = p % pRNG % get() ! Random number to sample sites ! Retrieve cross section at the energy used for reaction sampling - call self % nuc % getMicroXSs(microXSs, collDat % E, self % mat % kT, p % pRNG) + if (p % T <= ZERO) then + kT = self % mat % kT + else + kT = p % T * kBoltzmannMeV + end if + call self % nuc % getMicroXSs(microXSs, collDat % E, kT, p % pRNG) sig_nufiss = microXSs % nuFission sig_fiss = microXSs % fission @@ -561,6 +576,7 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) else collDat % kT = self % nuc % getkT() end if + if (p % T > ZERO) collDat % kT = p % T * kBoltzmannMeV ! Check is DBRC is on hasDBRC = self % nuc % hasDBRC() diff --git a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 index c2e4e757a..e1154bde2 100644 --- a/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 +++ b/CollisionOperator/CollisionProcessors/neutronCEstd_class.f90 @@ -2,7 +2,7 @@ module neutronCEstd_class use numPrecision use endfConstants - use universalVariables, only : REJECTED, INF + use universalVariables, only : REJECTED, kBoltzmannMeV use genericProcedures, only : fatalError, rotateVector, numToChar use dictionary_class, only : dictionary use RNG_class, only : RNG @@ -162,7 +162,7 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) class(particleDungeon),intent(inout) :: thisCycle class(particleDungeon),intent(inout) :: nextCycle type(neutronMicroXSs) :: microXSs - real(defReal) :: r, denom, alphaXS, probAlpha + real(defReal) :: r, kT, denom, alphaXS, probAlpha character(100),parameter :: Here = 'sampleCollision (neutronCEstd_class.f90)' ! Verify that particle is CE neutron @@ -191,10 +191,11 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) ! Verify and load material pointer self % mat => ceNeutronMaterial_CptrCast(self % xsData % getMaterial(p % matIdx())) - if (.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial') - + if (.not.associated(self % mat)) call fatalError(Here, 'Material is not ceNeutronMaterial: '& + //numToChar(p % matIdx())) + ! Select collision nuclide - call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E) + call self % mat % sampleNuclide(p % E, p % pRNG, collDat % nucIdx, collDat % E, p % T, p % rho) ! If nuclide was rejected in TMS loop return to tracking if (collDat % nucIdx == REJECTED) then @@ -206,7 +207,12 @@ subroutine sampleCollision(self, p, tally, collDat, thisCycle, nextCycle) if (.not.associated(self % nuc)) call fatalError(Here, 'Failed to retrieve CE Neutron Nuclide') ! Select Main reaction channel - call self % nuc % getMicroXSs(microXss, collDat % E, self % mat % kT, p % pRNG) + if (p % T <= ZERO) then + kT = self % mat % kT + else + kT = p % T * kBoltzmannMeV + end if + call self % nuc % getMicroXSs(microXss, collDat % E, kT, p % pRNG) r = p % pRNG % get() collDat % MT = microXss % invert(r) @@ -228,7 +234,7 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) real(defReal),dimension(3) :: r, dir integer(shortInt) :: n, i real(defReal) :: wgt, w0, rand1, E_out, mu, phi - real(defReal) :: sig_nufiss, sig_tot, k_eff, lambda, wD + real(defReal) :: sig_nufiss, sig_tot, k_eff, kT, lambda, wD character(100),parameter :: Here = 'implicit (neutronCEstd_class.f90)' ! Generate fission sites if nuclide is fissile @@ -241,7 +247,12 @@ subroutine implicit(self, p, tally, collDat, thisCycle, nextCycle) rand1 = p % pRNG % get() ! Random number to sample sites ! Retrieve cross section at the energy used for reaction sampling - call self % nuc % getMicroXSs(microXSs, collDat % E, self % mat % kT, p % pRNG) + if (p % T <= ZERO) then + kT = self % mat % kT + else + kT = p % T * kBoltzmannMeV + end if + call self % nuc % getMicroXSs(microXSs, collDat % E, kT, p % pRNG) sig_nufiss = microXSs % nuFission sig_tot = microXSs % total @@ -367,6 +378,7 @@ subroutine elastic(self, p, tally, collDat, thisCycle, nextCycle) else collDat % kT = self % nuc % getkT() end if + if (p % T > ZERO) collDat % kT = p % T * kBoltzmannMeV ! Check is DBRC is on hasDBRC = self % nuc % hasDBRC() diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 index c262e7ea8..e1b109513 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/cartesianField_class.f90 @@ -57,16 +57,16 @@ module cartesianField_class real(defReal), dimension(3) :: corner = ZERO real(defReal), dimension(3) :: a_bar = ZERO type(box) :: outline - integer(shortInt) :: outLocalID = 0 integer(shortInt) :: nMat = 0 integer(shortInt), dimension(:), allocatable :: matIdxs contains ! Superclass procedures - procedure :: init + procedure :: init_dict procedure :: at procedure :: atP procedure :: distance + procedure :: map procedure :: kill procedure, private :: getLocalID @@ -77,7 +77,7 @@ module cartesianField_class !! !! Initialisation !! - subroutine init(self, dict) + subroutine init_dict(self, dict) class(cartesianField), intent(inout) :: self class(dictionary), intent(in) :: dict type(dictionary) :: tempDict @@ -163,9 +163,8 @@ subroutine init(self, dict) end if ! Size field value array - self % outLocalID = product(self % sizeN) + 1 - self % N = product(self % sizeN * self % nMat) + 1 - allocate(self % val(self % N)) + self % N = product(self % sizeN * self % nMat) + allocate(self % val(self % N + 1)) ! Read field values for each material idx0 = 0 @@ -195,9 +194,9 @@ subroutine init(self, dict) end do ! Set default value when not in the field - call dict % getOrDefault(self % val(self % N), 'default', -INF) + call dict % getOrDefault(self % val(self % N + 1), 'default', -INF) - end subroutine init + end subroutine init_dict !! !! Get value of the field at the co-ordinate point @@ -208,25 +207,43 @@ function at(self, coords) result(val) class(cartesianField), intent(in) :: self class(coordList), intent(in) :: coords real(defReal) :: val - integer(shortInt) :: idx, localID + integer(shortInt) :: localID - localID = self % getLocalID(coords % lvl(1) % r, coords % lvl(1) % dir) + localID = self % map (coords) + if (localID == 0) localID = self % N + 1 + val = self % val(localID) - if (localID == self % outLocalID) then - val = self % val(self % N) + end function at + + !! + !! Get index of the field at the co-ordinate point + !! + !! See pieceConstantField for details + !! + pure function map(self, coords) result(idx) + class(cartesianField), intent(in) :: self + class(coordList), intent(in) :: coords + integer(shortInt) :: idx + integer(shortInt) :: idx0 + + idx = self % getLocalID(coords % lvl(1) % r, coords % lvl(1) % dir) + + ! Outside the field + if (idx == 0) then return end if ! Compare against material idx - if (self % matIdxs(1) /= ALL_MATS) then - idx = findloc(self % matIdxs, coords % matIdx, 1) - localID = localID + (idx - 1) * product(self % sizeN) + ! Ensure material is present + if (any(self % matIdxs == coords % matIdx)) then + idx0 = findloc(self % matIdxs, coords % matIdx, 1) + idx = idx + (idx0 - 1) * product(self % sizeN) + else if (self % matIdxs(1) /= ALL_MATS) then + idx = 0 end if - - val = self % val(localID) - end function at + end function map !! !! Get value of the field at the particle's location @@ -259,11 +276,17 @@ function distance(self, coords) result(d) localID = self % getLocalID(coords % lvl(1) % r, coords % lvl(1) % dir) ! Catch case if particle is outside the lattice - if (localID == self % outLocalID) then + if (localID == 0) then d = self % outline % distance(coords % lvl(1) % r, coords % lvl(1) % dir) return end if + + ! Catch case if particle is in an excluded material + if ((.not. any(coords % matIdx == self % matIdxs)) .and. self % matIdxs(1) /= ALL_MATS) then + d = INF + return + end if ! Compute ijk of localID temp = localID - 1 @@ -306,8 +329,7 @@ function distance(self, coords) result(d) end do ! Cap distance value - d = max(ZERO, d) - d = min(INF, d) + if (d <= ZERO) d = INF end function distance @@ -325,7 +347,6 @@ elemental subroutine kill(self) self % nMat = 0 self % corner = ZERO self % a_bar = ZERO - self % outLocalID = 0 call self % outline % kill() end subroutine kill @@ -333,14 +354,14 @@ end subroutine kill !! !! Find the local integer ID in the field given position and direction !! - function getLocalID(self, r, u) result(localID) - class(cartesianField), intent(in) :: self - real(defReal), dimension(3) :: r - real(defReal), dimension(3) :: u - integer(shortInt) :: localID - real(defReal), dimension(3) :: r_bar - integer(shortInt), dimension(3) :: ijk - integer(shortInt) :: i, inc + pure function getLocalID(self, r, u) result(localID) + class(cartesianField), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + integer(shortInt) :: localID + real(defReal), dimension(3) :: r_bar + integer(shortInt), dimension(3) :: ijk + integer(shortInt) :: i, inc ijk = floor((r - self % corner) / self % pitch) + 1 @@ -364,7 +385,7 @@ function getLocalID(self, r, u) result(localID) end do if (any(ijk <= 0 .or. ijk > self % sizeN)) then ! Point is outside lattice - localID = self % outLocalID + localID = 0 else localID = ijk(1) + self % sizeN(1) * (ijk(2)-1 + self % sizeN(2) * (ijk(3)-1)) diff --git a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 index 3f00690f5..8d21708e5 100644 --- a/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 +++ b/Geometry/Fields/ScalarFields/PieceConstantFields/pieceConstantField_inter.f90 @@ -5,6 +5,8 @@ module pieceConstantField_inter use scalarField_inter, only : scalarField use coord_class, only : coordList use dictionary_class, only : dictionary + use particle_class, only : particle + use dictParser_func, only : fileToDict use field_inter, only : field use scalarField_inter, only : scalarField @@ -22,26 +24,52 @@ module pieceConstantField_inter !! Provides a distance calculation to neighbouring elements, allowing surface !! tracking to be used. !! + !! Access to field is via coordList to allow more fancy fields to be defined + !! (e.g. assign value to each uniqueID etc.) + !! !! Interface: !! field interface - !! distance -> Return distance to next element of the field - !! setValue -> Sets value of field at a given index - !! getMaxValue -> Returns the maximum field value + !! init -> Initialise from dictionary, which may contain a path + !! init_dict -> Initialises from a dictionary containing full description + !! init_file -> Initialises from a path to a full description + !! map -> Returns field index given coordList + !! map_particle -> Returns field index given particle + !! distance -> Return distance to next element of the field + !! setValue -> Sets value of field at a given index + !! getMaxValue -> Returns the maximum field value + !! getSize -> Returns the size of the field !! type, public, abstract, extends(scalarField) :: pieceConstantField real(defReal), dimension(:), allocatable :: val integer(shortInt) :: N = 0 contains - procedure(distance), deferred :: distance + procedure(init_dict), deferred, private :: init_dict + procedure(distance), deferred :: distance + procedure(map), deferred :: map + procedure :: init + procedure :: map_particle procedure :: setValue procedure :: getMaxValue + procedure :: getSize procedure :: kill end type pieceConstantField abstract interface + !! + !! Initialise field from a dictionary + !! + !! Args: + !! dict [in] -> Dictionary describing field + !! + subroutine init_dict(self, dict) + import :: pieceConstantField, dictionary + class(pieceConstantField), intent(inout) :: self + class(dictionary), intent(in) :: dict + end subroutine init_dict + !! !! Get distance to the next element of the field at the given coordinates !! @@ -58,10 +86,49 @@ function distance(self, coords) result(d) real(defReal) :: d end function distance + !! + !! Get index of the field at the given coordinates + !! + !! Args: + !! coords [in] -> Coordinates of the position in the geometry + !! + !! Result: + !! Index into the field. Integer + !! + pure function map(self, coords) result(idx) + import :: pieceConstantField, coordList, shortInt + class(pieceConstantField), intent(in) :: self + class(coordList), intent(in) :: coords + integer(shortInt) :: idx + end function map + end interface contains + !! + !! Initialise field from a dictionary. This dictionary may contain + !! either the details of the field or a path to a further dictionary + !! + !! Args: + !! dict [in] -> Dictionary describing contents of field + !! + subroutine init(self, dict) + class(pieceConstantField), intent(inout) :: self + class(dictionary), intent(in) :: dict + character(pathLen) :: path + class(dictionary), pointer :: tempDict + + if (dict % isPresent('path')) then + call dict % get(path, 'path') + call fileToDict(tempDict, path) + call self % init_dict(tempDict) + else + call self % init_dict(dict) + end if + + end subroutine init + !! !! Cast field pointer to pieceConstantField pointer !! @@ -114,6 +181,20 @@ subroutine setValue(self, val, idx) self % val(idx) = val end subroutine setValue + + !! + !! Returns the index of the field given a particle. + !! + function map_particle(self, p) result(idx) + class(pieceConstantField), intent(in) :: self + type(particle), intent(in) :: p + integer(shortInt) :: idx + type(coordList) :: coords + + coords = p % coords + idx = self % map(coords) + + end function map_particle !! !! Returns the value of the maximum value of the field. @@ -123,9 +204,9 @@ end subroutine setValue !! If the values of the field are not allocated !! function getMaxValue(self) result(val) - class(pieceConstantField), intent(inout) :: self - real(defReal) :: val - character(100), parameter :: Here = 'getMaxValue (pieceConstantField_class.f90)' + class(pieceConstantField), intent(in) :: self + real(defReal) :: val + character(100), parameter :: Here = 'getMaxValue (pieceConstantField_class.f90)' if (.not. allocated(self % val)) call fatalError(Here,& 'The values of the field have not been allocated') @@ -135,6 +216,16 @@ function getMaxValue(self) result(val) end function getMaxValue !! + !! Returns the size of the field. + !! + function getSize(self) result(N) + class(pieceConstantField), intent(in) :: self + integer(shortInt) :: N + + N = self % N + + end function getSize + !! Kills elements of the superclass !! elemental subroutine kill(self) diff --git a/Geometry/Tests/geometryStd_iTest.f90 b/Geometry/Tests/geometryStd_iTest.f90 index 75def8e56..976e821ac 100644 --- a/Geometry/Tests/geometryStd_iTest.f90 +++ b/Geometry/Tests/geometryStd_iTest.f90 @@ -7,6 +7,7 @@ module geometryStd_iTest use dictParser_func, only : fileToDict use coord_class, only : coordList use geometryStd_class, only : geometryStd + use materialMenu_mod, only : mm_init => init use funit implicit none @@ -294,5 +295,4 @@ subroutine test_tilted_cylinder() end subroutine test_tilted_cylinder - end module geometryStd_iTest diff --git a/Geometry/coord_class.f90 b/Geometry/coord_class.f90 index b578254b4..fe8376f9d 100644 --- a/Geometry/coord_class.f90 +++ b/Geometry/coord_class.f90 @@ -300,7 +300,7 @@ elemental subroutine takeAboveGeom(self) end subroutine takeAboveGeom !! - !! Decrease nestting to level n + !! Decrease nesting to level n !! !! Args: !! n [in] -> New nesting level diff --git a/Geometry/fieldFactory_func.f90 b/Geometry/fieldFactory_func.f90 index 1d587072b..3e292e1c2 100644 --- a/Geometry/fieldFactory_func.f90 +++ b/Geometry/fieldFactory_func.f90 @@ -13,6 +13,7 @@ module fieldFactory_func use uniformVectorField_class, only : uniformVectorField use uniFissSitesField_class, only : uniFissSitesField use weightWindowsField_class, only : weightWindowsField + use cartesianField_class, only : cartesianField ! Geometry use geometryReg_mod, only : gr_addField => addField @@ -25,7 +26,8 @@ module fieldFactory_func character(nameLen), dimension(*), parameter :: AVAILABLE_FIELDS = ['uniformScalarField',& 'uniformVectorField',& 'uniFissSitesField ',& - 'weightWindowsField'] + 'weightWindowsField',& + 'cartesianField '] ! Public interface public :: new_field @@ -41,7 +43,7 @@ module fieldFactory_func !! name [in] -> Name of the field for the geometry registry !! !! Errors: - !! fatalError is type of field is unknown + !! fatalError if type of field is unknown !! subroutine new_field(dict, name) class(dictionary), intent(in) :: dict @@ -67,6 +69,9 @@ subroutine new_field(dict, name) case ('weightWindowsField') allocate(weightWindowsField :: kentta) + case ('cartesianField') + allocate(cartesianField :: kentta) + case default print '(A)', "AVAILABLE FIELDS:" print '(A)', AVAILABLE_FIELDS @@ -81,6 +86,5 @@ subroutine new_field(dict, name) call gr_addField(kentta, name) end subroutine new_field - - + end module fieldFactory_func diff --git a/Geometry/geometryReg_mod.f90 b/Geometry/geometryReg_mod.f90 index 153bd4443..a4ceb805b 100644 --- a/Geometry/geometryReg_mod.f90 +++ b/Geometry/geometryReg_mod.f90 @@ -16,27 +16,30 @@ !! fieldNameMap -> Map of field names to idx !! !! Interface: -!! addGeom -> Add new geometry -!! geomIdx -> Get index of the geometry -!! geomPtr -> Get pointer to a geometry with an index -!! addField -> Add new field -!! fieldIdx -> Get index of a field -!! fieldPtr -> Get pointer to a field specified by index -!! display -> Display info abou defined fields and geometries -!! kill -> Return to uninitialised state +!! addGeom -> Add new geometry +!! geomIdx -> Get index of the geometry +!! geomPtr -> Get pointer to a geometry with an index +!! addField -> Add new field +!! hasField -> Does the field with a given name exist? +!! fieldIdx -> Get index of a field +!! fieldPtr -> Get pointer to a field specified by index +!! fieldPtrName -> Get pointer to a field specified by name, for select names +!! display -> Display info about defined fields and geometries +!! kill -> Return to uninitialised state !! module geometryReg_mod use numPrecision - use genericProcedures, only : fatalError, numToChar - use dictionary_class, only : dictionary - use charMap_class, only : charMap + use genericProcedures, only : fatalError, numToChar + use universalVariables, only : nameTemperature, nameDensity + use dictionary_class, only : dictionary + use charMap_class, only : charMap ! Geometry - use geometry_inter, only : geometry + use geometry_inter, only : geometry ! Fields - use field_inter, only : field + use field_inter, only : field implicit none private @@ -66,6 +69,8 @@ module geometryReg_mod public :: addField public :: fieldIdx public :: fieldPtr + public :: fieldPtrName + public :: hasField public :: kill integer(shortInt), parameter :: START_SIZE = 5 @@ -79,6 +84,10 @@ module geometryReg_mod type(fieldBox), dimension(:), allocatable, target :: fields integer(shortInt) :: fieldTop = 0 type(charMap) :: fieldNameMap + + ! Save indices for speed. Avoids repeated calls to the character map + integer(shortInt) :: temperatureIdx = -8 + integer(shortInt) :: densityIdx = -8 contains @@ -226,10 +235,45 @@ subroutine addField(kentta, name) call fieldNameMap % add(name, idx) fields(idx) % name = name + ! Save density or temperature indices for speed + if (name == nameTemperature) then + temperatureIdx = idx + elseif (name == nameDensity) then + densityIdx = idx + end if + ! Point field call move_alloc(kentta, fields(idx) % kentta) end subroutine addField + + !! + !! Returns whether a field is present + !! + !! Args: + !! name [in] -> Name of the field + !! + !! Result: + !! Logical stating whether the field exists + !! + pure function hasField(name) result(exists) + character(nameLen), intent(in) :: name + logical(defBool) :: exists + integer(shortInt), parameter :: NOT_PRESENT = -8 + integer(shortInt) :: idx + + select case(name) + case(nameTemperature) + idx = temperatureIdx + case(nameDensity) + idx = densityIdx + case default + idx = fieldNameMap % getOrDefault(name, NOT_PRESENT) + end select + + exists = (idx /= NOT_PRESENT) + + end function hasField !! !! Get index of a field given its name @@ -255,6 +299,46 @@ function fieldIdx(name) result(idx) end if end function fieldIdx + + !! + !! Get pointer to a field given its name + !! Optimised for frequently accessed fields to avoid + !! the relatively slow character map + !! + !! Args: + !! name [in] -> Name of the field + !! + !! Result: + !! Pointer to the field with the name + !! + !! Errors: + !! fatalError if name was not defined or + !! an unrecognised name was given + !! + function fieldPtrName(name) result(ptr) + character(nameLen), intent(in) :: name + class(field), pointer :: ptr + integer(shortInt) :: idx + integer(shortInt), parameter :: NOT_PRESENT = -8 + character(*), parameter :: Here = 'fieldPtrName (geometryReg_mod.f90)' + + select case(name) + case(nameTemperature) + idx = temperatureIdx + case(nameDensity) + idx = densityIdx + case default + idx = NOT_PRESENT + end select + + if (idx < 1 .or. idx > fieldTop) then + call fatalError(Here,'Index: '//numToChar(idx)//' does not correspond to valid field. & + &Must be 1-'//numToChar(fieldTop)) + end if + + ptr => fields(idx) % kentta + + end function fieldPtrName !! !! Get pointer to a field by index diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index 49fb71230..d3a0e03f4 100644 --- a/Geometry/geometryStd_class.f90 +++ b/Geometry/geometryStd_class.f90 @@ -10,6 +10,12 @@ module geometryStd_class use csg_class, only : csg use universe_inter, only : universe use surface_inter, only : surface + + use field_inter, only : field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast + + use geometryReg_mod, only : gr_hasField => hasField, & + gr_fieldPtrName => fieldPtrName ! Nuclear Data use materialMenu_mod, only : nMat @@ -23,12 +29,19 @@ module geometryStd_class !! public :: geometryStd_CptrCast + !! + !! Private distance to the nearest pieceConstant field crossing + !! + private :: getFieldDist + !! !! Standard Geometry Model !! !! Typical geometry of a MC Neutron Transport code composed of multiple nested !! universes. !! + !! Supports super-imposed temperature and density fields. + !! !! Boundary conditions in diffrent movement models are handeled: !! move -> explicitBC !! moveGlobal -> explicitBC @@ -37,7 +50,7 @@ module geometryStd_class !! Sample Dictionary Input: !! geometry { !! type geometryStd; - !! + !! !! } !! !! Public Members: @@ -198,28 +211,38 @@ end function bounds !! Uses explicit BC !! subroutine move_noCache(self, coords, maxDist, event) - class(geometryStd), intent(in) :: self - type(coordList), intent(inout) :: coords - real(defReal), intent(inout) :: maxDist - integer(shortInt), intent(out) :: event - integer(shortInt) :: surfIdx, level - real(defReal) :: dist - class(surface), pointer :: surf - class(universe), pointer :: uni - character(100), parameter :: Here = 'move (geometryStd_class.f90)' + class(geometryStd), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + integer(shortInt) :: surfIdx, level, level0 + real(defReal) :: dist, fieldDist + class(surface), pointer :: surf + class(universe), pointer :: uni + character(100), parameter :: Here = 'move_noCache (geometryStd_class.f90)' if (.not.coords % isPlaced()) then call fatalError(Here, 'Coordinate list is not placed in the geometry') end if + level0 = coords % nesting + ! Find distance to the next surface call self % closestDist(dist, surfIdx, level, coords) + + ! Check fields + fieldDist = getFieldDist(coords) - if (maxDist < dist) then ! Moves within cell + if (maxDist < dist .and. maxDist < fieldDist) then ! Moves within cell call coords % moveLocal(maxDist, coords % nesting) event = COLL_EV maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway - + + else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + call coords % moveLocal(fieldDist, level0) + event = FIELD_EV + maxDist = fieldDist + else if (surfIdx == self % geom % borderIdx .and. level == 1) then ! Hits domain boundary ! Move global to the boundary call coords % moveGlobal(dist) @@ -232,7 +255,7 @@ subroutine move_noCache(self, coords, maxDist, event) ! Place back in geometry call self % placeCoord(coords) - + else ! Crosses to different local cell ! Move to boundary at hit level call coords % moveLocal(dist, level) @@ -248,6 +271,7 @@ subroutine move_noCache(self, coords, maxDist, event) end if + end subroutine move_noCache !! @@ -258,30 +282,41 @@ end subroutine move_noCache !! Uses explicit BC !! subroutine move_withCache(self, coords, maxDist, event, cache) - class(geometryStd), intent(in) :: self - type(coordList), intent(inout) :: coords - real(defReal), intent(inout) :: maxDist - integer(shortInt), intent(out) :: event - type(distCache), intent(inout) :: cache - integer(shortInt) :: surfIdx, level - real(defReal) :: dist - class(surface), pointer :: surf - class(universe), pointer :: uni + class(geometryStd), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + type(distCache), intent(inout) :: cache + integer(shortInt) :: surfIdx, level, level0 + real(defReal) :: dist, fieldDist + class(surface), pointer :: surf + class(universe), pointer :: uni character(100), parameter :: Here = 'move_withCache (geometryStd_class.f90)' if (.not.coords % isPlaced()) then call fatalError(Here, 'Coordinate list is not placed in the geometry') end if + level0 = coords % nesting + ! Find distance to the next surface call self % closestDist_cache(dist, surfIdx, level, coords, cache) + + ! Check fields + fieldDist = getFieldDist(coords) - if (maxDist < dist) then ! Moves within cell + if (maxDist < dist .and. maxDist < fieldDist) then ! Moves within cell call coords % moveLocal(maxDist, coords % nesting) event = COLL_EV maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway cache % lvl = 0 + else if (maxDist < dist .and. maxDist >= fieldDist) then ! Stays within the same cell, but crosses field boundary + call coords % moveLocal(fieldDist, level0) + event = FIELD_EV + maxDist = fieldDist + cache % dist(1:level0) = cache % dist(1:level0) - fieldDist + else if (surfIdx == self % geom % borderIdx .and. level == 1) then ! Hits domain boundary ! Move global to the boundary call coords % moveGlobal(dist) @@ -295,7 +330,7 @@ subroutine move_withCache(self, coords, maxDist, event, cache) ! Place back in geometry call self % placeCoord(coords) - + else ! Crosses to different local cell ! Move to boundary at hit level call coords % moveLocal(dist, level) @@ -606,5 +641,33 @@ pure function geometryStd_CptrCast(source) result(ptr) end select end function geometryStd_CptrCast + + !! + !! Helper function. + !! Get the distance to the nearest pieceConstantField boundary + !! + !! Returns INF if neither density nor temperature fields exist + !! + function getFieldDist(coords) result(dist) + type(coordList), intent(in) :: coords + real(defReal) :: dist + class(field), pointer :: genericField + class(pieceConstantField), pointer :: pcField + + dist = INF + if (gr_hasField(nameTemperature)) then + genericField => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(genericField) + dist = min(dist, pcField % distance(coords)) + end if + + if (gr_hasField(nameDensity)) then + genericField => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(genericField) + dist = min(dist, pcField % distance(coords)) + end if + + end function getFieldDist + end module geometryStd_class diff --git a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 index cb2001baa..d5717486d 100644 --- a/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/Tests/aceNeutronDatabase_iTest.f90 @@ -237,7 +237,7 @@ subroutine test_aceNeutronDatabase() ! ! Water mat => ceNeutronMaterial_TptrCast( data % getMaterial(1)) - call mat % getMacroXSs(macroXss, 3.6E-1_defReal, p % pRNG) + call mat % getMacroXSs(macroXss, 3.6E-1_defReal, p % T, p % rho, p % pRNG) ! Absent XSs @assertEqual(ZERO, macroXSs % fission) @@ -249,7 +249,7 @@ subroutine test_aceNeutronDatabase() @assertEqual(ONE, 2.198066842597500e-06_defReal/ macroXSs % capture, TOL) ! Water with some inelastic collisions - call mat % getMacroXSs(macroXss, 6.525_defReal, p % pRNG) + call mat % getMacroXSs(macroXss, 6.525_defReal, p % T, p % rho, p % pRNG) @assertEqual(ONE, macroXSs % inelasticScatter/1.903667536E-04_defReal, TOL) diff --git a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 index 6899dd076..e503fc6a5 100644 --- a/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 +++ b/NuclearData/ceNeutronData/aceDatabase/aceNeutronDatabase_class.f90 @@ -103,6 +103,7 @@ module aceNeutronDatabase_class procedure :: getReaction procedure :: init procedure :: activate + procedure :: initMajorant ! ceNeutronDatabase Procedures procedure :: energyBounds @@ -119,7 +120,6 @@ module aceNeutronDatabase_class ! class Procedures procedure :: initUrr procedure :: initDBRC - procedure :: initMajorant procedure :: updateTotalTempMajXS procedure :: updateRelEnMacroXSs procedure :: makeNuclideName @@ -370,7 +370,8 @@ subroutine updateMajorantXS(self, E, rand) maj % xs = self % majorant(idx+1) * f + (ONE - f) * self % majorant(idx) - else ! Compute majorant on the fly + ! Compute majorant on the fly - will fail if using a density map! + else maj % xs = ZERO @@ -380,7 +381,7 @@ subroutine updateMajorantXS(self, E, rand) ! Update if needed if (cache_materialCache(matIdx) % E_track /= E) then - call self % updateTrackMatXS(E, matIdx, rand) + call self % updateTrackMatXS(E, matIdx, rand = rand) end if maj % xs = max(maj % xs, cache_materialCache(matIdx) % trackXS) @@ -398,25 +399,41 @@ end subroutine updateMajorantXS !! !! See ceNeutronDatabase for more details !! - subroutine updateTrackMatXS(self, E, matIdx, rand) + subroutine updateTrackMatXS(self, E, matIdx, temp, rho, rand) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), optional, intent(in) :: temp + real(defReal), optional, intent(in) :: rho class(RNG), optional, intent(inout) :: rand + real(defReal) :: T, densityFactor associate (matCache => cache_materialCache(matIdx), & mat => self % materials(matIdx)) - ! Set new energy + ! Set new energy, temperature, and density factor matCache % E_track = E + if (present(temp)) then + T = temp + else + T = -ONE + end if + matCache % T_track = T + + if (present(rho)) then + densityFactor = rho + else + densityFactor = -ONE + end if + matCache % rho_track = densityFactor if (mat % useTMS(E)) then ! The material tracking xs is the temperature majorant in the case of TMS - call self % updateTotalTempMajXS(E, matIdx) + call self % updateTotalTempMajXS(E, matIdx, T, densityFactor) else ! When TMS is not in use, the material tracking xs is equivalent to the total - call self % updateTotalMatXS(E, matIdx, rand) + call self % updateTotalMatXS(E, matIdx, T, densityFactor, rand) matCache % trackXS = matCache % xss % total end if @@ -436,12 +453,14 @@ end subroutine updateTrackMatXS !! E [in] -> Incident neutron energy for which temperature majorant is found !! matIdx [in] -> Index of material for which the material temperature majorant is found !! - subroutine updateTotalTempMajXS(self, E, matIdx) + subroutine updateTotalTempMajXS(self, E, matIdx, temp, rho) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho integer(shortInt) :: nucIdx, i - real(defReal) :: dens, corrFact, nucTempMaj + real(defReal) :: kT, dens, corrFact, nucTempMaj associate (matCache => cache_materialCache(matIdx), & mat => self % materials(matIdx)) @@ -449,6 +468,13 @@ subroutine updateTotalTempMajXS(self, E, matIdx) ! Clean current total XS matCache % trackXS = ZERO + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = mat % kT + else + kT = temp * kBoltzmannMeV + end if + ! loop through all nuclides in material and find sum of majorants do i = 1, size(mat % nuclides) @@ -456,7 +482,7 @@ subroutine updateTotalTempMajXS(self, E, matIdx) nucIdx = mat % nuclides(i) dens = mat % dens(i) - call self % updateTotalTempNucXS(E, mat % kT, nucIdx) + call self % updateTotalTempNucXS(E, kT, nucIdx) ! Sum nuclide majorants to find material majorant corrFact = cache_nuclideCache(nucIdx) % doppCorr @@ -465,6 +491,11 @@ subroutine updateTotalTempMajXS(self, E, matIdx) end do + ! Use imposed density scaling if given + if (rho > ZERO) then + matCache % trackXS = matCache % trackXS * rho + end if + end associate end subroutine updateTotalTempMajXS @@ -475,26 +506,43 @@ end subroutine updateTotalTempMajXS !! !! See ceNeutronDatabase for more details !! - subroutine updateTotalMatXS(self, E, matIdx, rand) + subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), optional, intent(inout) :: rand integer(shortInt) :: i, nucIdx - real(defReal) :: dens + real(defReal) :: dens, kT, densityFactor associate (matCache => cache_materialCache(matIdx), & mat => self % materials(matIdx)) ! Set new energy and clean current total XS - matCache % E_tot = E matCache % xss % total = ZERO + matCache % E_tot = E + matCache % rho_tot = rho + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = mat % kT + else + kT = temp * kBoltzmannMeV + end if + + ! Use imposed density scaling if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if if (mat % useTMS(E)) then ! When TMS is in use, the total xs is retrieved sampling the nuclides' relative ! energies given the temperature difference between material temperature and ! temperature of the nuclides' base cross sections - call self % updateRelEnMacroXSs(E, matIdx, rand) + call self % updateRelEnMacroXSs(E, matIdx, kT, densityFactor, rand) else ! Construct total macro XS @@ -503,16 +551,20 @@ subroutine updateTotalMatXS(self, E, matIdx, rand) nucIdx = mat % nuclides(i) ! Update if needed - if (cache_nuclideCache(nucIdx) % E_tot /= E) then - call self % updateTotalNucXS(E, nucIdx, mat % kT, rand) + if (cache_nuclideCache(nucIdx) % E_tot /= E .or. matCache % T_tot /= temp) then + call self % updateTotalNucXS(E, nucIdx, kT, rand) end if ! Add microscopic XSs matCache % xss % total = matCache % xss % total + & dens * cache_nuclideCache(nucIdx) % xss % total end do + + matCache % xss % total = matCache % xss % total * densityFactor end if + + matCache % T_tot = temp end associate @@ -524,25 +576,42 @@ end subroutine updateTotalMatXS !! !! See ceNeutronDatabase for more details !! - subroutine updateMacroXSs(self, E, matIdx, rand) + subroutine updateMacroXSs(self, E, matIdx, temp, rho, rand) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), optional, intent(inout) :: rand integer(shortInt) :: i, nucIdx - real(defReal) :: dens + real(defReal) :: dens, kT, densityFactor associate(mat => self % materials(matIdx), & matCache => cache_materialCache(matIdx)) ! Clean current xss call matCache % xss % clean() + matCache % rho_tot = rho + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = mat % kT + else + kT = temp * kBoltzmannMeV + end if + + ! Use imposed density scaling if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if if (mat % useTMS(E)) then ! When TMS is in use, the xss are retrieved sampling the nuclides' relative ! energies given the temperature difference between material temperature and ! temperature of the nuclides' base cross sections - call self % updateRelEnMacroXSs(E, matIdx, rand) + call self % updateRelEnMacroXSs(E, matIdx, kT, densityFactor, rand) else @@ -556,15 +625,19 @@ subroutine updateMacroXSs(self, E, matIdx, rand) nucIdx = mat % nuclides(i) ! Update if needed - if (cache_nuclideCache(nucIdx) % E_tail /= E .or. cache_nuclideCache(nucIdx) % E_tot /= E) then - call self % updateMicroXSs(E, nucIdx, mat % kT, rand) + if (cache_nuclideCache(nucIdx) % E_tail /= E .or. & + cache_nuclideCache(nucIdx) % E_tot /= E .or. & + matCache % T_tot /= temp) then + call self % updateMicroXSs(E, nucIdx, kT, rand) end if ! Add microscopic XSs - call matCache % xss % add(cache_nuclideCache(nucIdx) % xss, dens) + call matCache % xss % add(cache_nuclideCache(nucIdx) % xss, dens * densityFactor) end do end if + + matCache % T_tot = temp end associate @@ -579,10 +652,12 @@ end subroutine updateMacroXSs !! E [in] -> Incident neutron energy for which the relative energy xss are found !! matIdx [in] -> Index of material for which the relative energy xss are found !! - subroutine updateRelEnMacroXSs(self, E, matIdx, rand) + subroutine updateRelEnMacroXSs(self, E, matIdx, kT, densityFactor, rand) class(aceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: kT + real(defReal), intent(in) :: densityFactor class(RNG), optional, intent(inout) :: rand integer(shortInt) :: i, nucIdx real(defReal) :: dens, nuckT, A, deltakT, eRel, eMin, & @@ -606,7 +681,11 @@ subroutine updateRelEnMacroXSs(self, E, matIdx, rand) nucIdx = mat % nuclides(i) nuckT = self % nuclides(nucIdx) % getkT() A = self % nuclides(nucIdx) % getMass() - deltakT = mat % kT - nuckT + deltakT = kT - nuckT + + ! Catch negative dkT - perhaps due to input field values + if (deltakT < ZERO) call fatalError(Here,& + 'Invalid input temperature: '//numToChar(kT/kBoltzmannMeV)) eRel = relativeEnergy_constXS(E, A, deltakT, rand) @@ -624,11 +703,11 @@ subroutine updateRelEnMacroXSs(self, E, matIdx, rand) ! Update if needed if (nucCache % E_tail /= eRel .or. nucCache % E_tot /= eRel) then - call self % updateMicroXSs(eRel, nucIdx, mat % kT, rand) + call self % updateMicroXSs(eRel, nucIdx, kT, rand) end if ! Add microscopic XSs - call matCache % xssRel % add(nucCache % xss, dens * doppCorr) + call matCache % xssRel % add(nucCache % xss, dens * doppCorr * densityFactor) end associate @@ -751,6 +830,10 @@ subroutine updateTotalTempNucXS(self, E, kT, nucIdx) A = nuc % getMass() deltakT = kT - nuckT + ! Catch negative dkT - perhaps due to input field values + if (deltakT < ZERO) call fatalError(Here,& + 'Invalid input temperatureL '//numToChar(kT/kBoltzmannMeV)) + ! Check if an update is required if (nucCache % E_maj /= E .or. nucCache % deltakT /= deltakT) then @@ -988,12 +1071,12 @@ subroutine init(self, dict, ptr, silent ) ! is bounded by Sab temperatures if (mat % nuclides(j) % sabMix) then sabT = self % nuclides(nucIdxs(j)) % getSabTBounds() - kT = mat % T * kBoltzmann / joulesPerMeV + kT = mat % T * kBoltzmannMev if ((kT < sabT(1)) .or. (kT > sabT(2))) call fatalError(Here,& 'Material temperature must be bounded by the provided S(alpha,beta) data. '//& - 'The material temperature is '//numToChar(kT * joulesPerMeV / kBoltzmann)//& - 'K while the data bounds are '//numToChar(sabT(1) * joulesPerMeV / kBoltzmann)//& - 'K and '//numToChar(sabT(2) * joulesPerMeV / kBoltzmann)//'K.') + 'The material temperature is '//numToChar(kT / kBoltzmannMeV)//& + 'K while the data bounds are '//numToChar(sabT(1) / kBoltzmannMeV)//& + 'K and '//numToChar(sabT(2) / kBoltzmannMeV)//'K.') end if end do @@ -1241,20 +1324,51 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, loud) + subroutine initMajorant(self, loud, maxTemp, scaleDensity) class(aceNeutronDatabase), intent(inout) :: self - logical(defBool), intent(in) :: loud + logical(defBool), intent(in), optional :: loud + real(defReal), intent(in), optional :: maxTemp + real(defReal), intent(in), optional :: scaleDensity + logical(defBool) :: isLoud real(defReal), dimension(:), allocatable :: tmpGrid integer(shortInt) :: i, j, k, matIdx, nNuc, nucIdx, isDone, & sizeGrid, eIdx, nucIdxLast, eIdxLast, & urrIdx type(intMap) :: nucSet real(defReal) :: eRef, eNuc, E, maj, trackXS, dens, urrMaj, & - nucXS, f, eMax, eMin + nucXS, f, eMax, eMin, densityFactor, broadenTemp class(RNG), allocatable :: rand integer(shortInt), parameter :: IN_SET = 1, NOT_PRESENT = 0 real(defReal), parameter :: NUDGE = 1.0e-06_defReal + if (present(loud)) then + isLoud = loud + else + isLoud = .false. + end if + + ! Check maxTemp is present and sensible + if (present(maxTemp)) then + if (maxTemp < ZERO) then + broadenTemp = -ONE + else + broadenTemp = maxTemp + end if + else + broadenTemp = -ONE + end if + + ! Check scaleDensity is present and sensible + if (present(scaleDensity)) then + if (scaleDensity < ONE) then + densityFactor = ONE + else + densityFactor = scaleDensity + end if + else + densityFactor = ONE + end if + ! Find the size of the unionised energy grid (with duplicates) ! Initialise size sizeGrid = 0 @@ -1422,6 +1536,7 @@ subroutine initMajorant(self, loud) end if ! Allocate unionised majorant + if (allocated(self % majorant)) deallocate(self % majorant) allocate(self % majorant(size(self % eGridUnion))) ! Initialise RNG needed to call update XS routines. The initial seed doesn't @@ -1450,12 +1565,16 @@ subroutine initMajorant(self, loud) matIdx = self % activeMat(j) ! Get material tracking cross section - call self % updateTrackMatXS(E, matIdx, rand) + if (self % materials(matIdx) % hasTMS .and. broadenTemp > ZERO) then + call self % updateTrackMatXS(E, matIdx, temp = broadenTemp, rho = densityFactor, rand = rand) + else + call self % updateTrackMatXS(E, matIdx, rho = densityFactor, rand = rand) + end if trackXS = cache_materialCache(matIdx) % trackXS ! Loop over nuclides to check and correct for ures do k = 1, size(self % materials(matIdx) % nuclides) - dens = self % materials(matIdx) % dens(k) + dens = self % materials(matIdx) % dens(k) * densityFactor nucIdx = self % materials(matIdx) % nuclides(k) associate (nuc => self % nuclides(nucIdx)) diff --git a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 index 935df7300..60c5402c8 100644 --- a/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 +++ b/NuclearData/ceNeutronData/ceNeutronCache_mod.f90 @@ -29,26 +29,34 @@ module ceNeutronCache_mod !! Structure that contains cached data for each CE Neutron Material !! !! Public Members: - !! E_tot -> Energy of the total XS in xss - !! E_tail -> Energy of all XSs in xss except total - !! f -> Interpolation factor for the nuclide at energy E_tot - !! idx -> Index on a nuclide grid for energy E_tot - !! xss -> Cached cross-section values - !! E_track -> Energy of the tracking xs - !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS - !! E_rel -> Base energy for which relative energy cross sections are found (for TMS) - !! xssRel -> Cached effective cross-section values at energy relative to E_rel (for TMS) + !! E_tot -> Energy of the total XS in xss + !! E_tail -> Energy of all XSs in xss except total + !! T_tot -> Temperature at which total XS was evaluated + !! rho_tot -> Density scaling at which total XS was evaluated + !! f -> Interpolation factor for the nuclide at energy E_tot + !! idx -> Index on a nuclide grid for energy E_tot + !! xss -> Cached cross-section values + !! E_track -> Energy of the tracking xs + !! T_track -> Temperature at which tracking XS was evaluated + !! rho_track -> Density scaling at which tracking XS was evaluated + !! trackXS -> Cached tracking xs; this can be different to xss % total when using TMS + !! E_rel -> Base energy for which relative energy cross sections are found (for TMS) + !! xssRel -> Cached effective cross-section values at energy relative to E_rel (for TMS) !! type, public :: cacheMatDat real(defReal) :: E_tot = ZERO + real(defReal) :: T_tot = -ONE + real(defReal) :: rho_tot = -ONE real(defReal) :: E_tail = ZERO real(defReal) :: f = ZERO integer(shortInt) :: idx = 0 type(neutronMacroXSs) :: xss ! Tracking data - real(defReal) :: E_track = ZERO - real(defReal) :: trackXS = ZERO + real(defReal) :: E_track = ZERO + real(defReal) :: trackXS = ZERO + real(defReal) :: T_track = -ONE + real(defReal) :: rho_track = -ONE ! TMS data real(defReal) :: E_rel = ZERO @@ -93,12 +101,12 @@ module ceNeutronCache_mod !! Structure that contains a cross section value !! !! Public Members: - !! E -> energy of the cross section - !! xs -> value of the cross section + !! E -> energy of the cross section + !! xs -> value of the cross section !! type, public :: cacheSingleXS - real(defReal) :: E = ZERO - real(defReal) :: xs = ZERO + real(defReal) :: E = ZERO + real(defReal) :: xs = ZERO end type cacheSingleXS !! diff --git a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 index 8ca9d1665..d9352adb9 100644 --- a/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 +++ b/NuclearData/ceNeutronData/ceNeutronDatabase_inter.f90 @@ -64,10 +64,10 @@ module ceNeutronDatabase_inter ! Procedures implemented by a specific CE Neutron Database procedure(updateMajorantXS), deferred :: updateMajorantXS - procedure(updateTotalMatXS), deferred :: updateTrackMatXS + procedure(updateTrackMatXS), deferred :: updateTrackMatXS procedure(updateTotalMatXS), deferred :: updateTotalMatXS procedure(updateMacroXSs), deferred :: updateMacroXSs - procedure(updateTotalXS), deferred :: updateTotalNucXS + procedure(updateTotalNucXS), deferred :: updateTotalNucXS procedure(updateMicroXSs), deferred :: updateMicroXSs procedure(updateTotalTempNucXS), deferred :: updateTotalTempNucXS procedure(energyBounds), deferred :: energyBounds @@ -99,9 +99,9 @@ end subroutine energyBounds !! !! Make sure that trackXS of material with matIdx is at energy E = E_track - !! in ceNeutronChache + !! in ceNeutronCache !! - !! The tracking xs correspons to the material total cross section unless TMS + !! The tracking xs corresponds to the material total cross section unless TMS !! is used. In that case, this is the material temperature majorant xs. !! !! Assume that call to this procedure implies that data is NOT up-to-date @@ -109,19 +109,23 @@ end subroutine energyBounds !! Args: !! E [in] -> required energy [MeV] !! matIdx [in] -> material index that needs to be updated + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! - subroutine updateTrackMatXS(self, E, matIdx, rand) + subroutine updateTrackMatXS(self, E, matIdx, temp, rho, rand) import :: ceNeutronDatabase, defReal, shortInt, RNG class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in), optional :: temp + real(defReal), intent(in), optional :: rho class(RNG), optional, intent(inout) :: rand end subroutine updateTrackMatXS !! !! Make sure that totalXS of material with matIdx is at energy E - !! in ceNeutronChache + !! in ceNeutronCache !! !! ANY CHANGE in ceNeutronChache is POSSIBLE !! E.G. All material XSs may be updated to energy E @@ -131,13 +135,17 @@ end subroutine updateTrackMatXS !! Args: !! E [in] -> required energy [MeV] !! matIdx [in] -> material index that needs to be updated + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! - subroutine updateTotalMatXS(self, E, matIdx, rand) + subroutine updateTotalMatXS(self, E, matIdx, temp, rho, rand) import :: ceNeutronDatabase, defReal, shortInt, RNG class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), optional, intent(inout) :: rand end subroutine updateTotalMatXS @@ -145,7 +153,7 @@ end subroutine updateTotalMatXS !! Make sure that the majorant of ALL Active materials is at energy E !! in ceNeutronChache !! - !! ANY CHANGE in ceNeutronChache is POSSIBLE + !! ANY CHANGE in ceNeutronCache is POSSIBLE !! E.G. All material XSs may be updated to energy E !! !! Assume that call to this procedure implies that data is NOT up-to-date @@ -165,7 +173,7 @@ end subroutine updateMajorantXS !! Make sure that the macroscopic XSs for the material with matIdx are set !! to energy E in ceNeutronCache !! - !! ANY CHANGE in ceNeutronChache is POSSIBLE + !! ANY CHANGE in ceNeutronCache is POSSIBLE !! E.G. Extra materials may be set to energy E as well !! !! Assume that call to this procedure implies that data is NOT up-to-date @@ -173,13 +181,17 @@ end subroutine updateMajorantXS !! Args: !! E [in] -> required energy [MeV] !! matIdx [in] -> material index that needs to be updated + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! - subroutine updateMacroXSs(self, E, matIdx, rand) + subroutine updateMacroXSs(self, E, matIdx, temp, rho, rand) import :: ceNeutronDatabase, defReal, shortInt, RNG class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: matIdx + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), optional, intent(inout) :: rand end subroutine updateMacroXSs @@ -187,8 +199,8 @@ end subroutine updateMacroXSs !! Make sure that totalXS of nuclide with nucIdx is at energy E !! in ceNeutronChache !! - !! ANY CHANGE in ceNeutronChache is POSSIBLE - !! E.G. All nuclid XSs may be updated to energy E + !! ANY CHANGE in ceNeutronCache is POSSIBLE + !! E.G. All nuclide XSs may be updated to energy E !! !! Assume that call to this procedure implies that data is NOT up-to-date !! @@ -198,14 +210,14 @@ end subroutine updateMacroXSs !! kT [in] -> thermal energy of material [MeV] !! rand [inout] -> random number generator !! - subroutine updateTotalXS(self, E, nucIdx, kT, rand) + subroutine updateTotalNucXS(self, E, nucIdx, kT, rand) import :: ceNeutronDatabase, defReal, shortInt, RNG class(ceNeutronDatabase), intent(in) :: self real(defReal), intent(in) :: E integer(shortInt), intent(in) :: nucIdx real(defReal), intent(in) :: kT class(RNG), optional, intent(inout) :: rand - end subroutine updateTotalXS + end subroutine updateTotalNucXS !! !! Make sure that the microscopic XSs for the nuclide with nucIdx are set @@ -345,8 +357,8 @@ function getTrackingXS(self, p, matIdx, what) result(xs) xs = max(xs, self % collisionXS) ! Update Cache - trackingCache(1) % E = p % E - trackingCache(1) % xs = xs + trackingCache(1) % E = p % E + trackingCache(1) % xs = xs end function getTrackingXS @@ -370,7 +382,7 @@ function getTrackMatXS(self, p, matIdx) result(xs) ! Check dynamic type of the particle if (p % isMG .or. p % type /= P_NEUTRON) then - call fatalError(Here, 'Dynamic type of the partcle is not CE Neutron but:'//p % typeToChar()) + call fatalError(Here, 'Dynamic type of the particle is not CE Neutron but:'//p % typeToChar()) end if ! Check that matIdx exists @@ -384,7 +396,11 @@ function getTrackMatXS(self, p, matIdx) result(xs) end if ! Check Cache and update if needed - if (materialCache(matIdx) % E_track /= p % E) call self % updateTrackMatXS(p % E, matIdx, p % pRNG) + if (materialCache(matIdx) % E_track /= p % E .or. & + materialCache(matIdx) % T_track /= p % T .or. & + materialCache(matIdx) % rho_track /= p % rho) then + call self % updateTrackMatXS(p % E, matIdx, p % T, p % rho, p % pRNG) + end if ! Return Cross-Section xs = materialCache(matIdx) % trackXS @@ -411,7 +427,7 @@ function getTotalMatXS(self, p, matIdx) result(xs) ! Check dynamic type of the particle if (p % isMG .or. p % type /= P_NEUTRON) then - call fatalError(Here, 'Dynamic type of the partcle is not CE Neutron but:'//p % typeToChar()) + call fatalError(Here, 'Dynamic type of the particle is not CE Neutron but:'//p % typeToChar()) end if ! Check that matIdx exists @@ -425,7 +441,11 @@ function getTotalMatXS(self, p, matIdx) result(xs) end if ! Check Cache and update if needed - if (materialCache(matIdx) % E_tot /= p % E) call self % updateTotalMatXS(p % E, matIdx, p % pRNG) + if (materialCache(matIdx) % E_tot /= p % E .or. & + materialCache(matIdx) % T_tot /= p % T .or. & + materialCache(matIdx) % rho_tot /= p % rho) then + call self % updateTotalMatXS(p % E, matIdx, p % T, p % rho, p % pRNG) + end if ! Return Cross-Section xs = materialCache(matIdx) % xss % total diff --git a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 index 0480c7034..05395c74f 100644 --- a/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 +++ b/NuclearData/ceNeutronData/ceNeutronMaterial_class.f90 @@ -115,8 +115,8 @@ subroutine getMacroXSs_byP(self, xss, p) class(particle), intent(in) :: p character(100), parameter :: Here = 'getMacroXSs_byP (ceNeutronMaterial_class.f90)' - if (.not. p % isMG) then - call self % getMacroXSs(xss, p % E, p % pRNG) + if (.not.p % isMG) then + call self % getMacroXSs(xss, p % E, p % T, p % rho, p % pRNG) else call fatalError(Here,'MG neutron given to CE data') @@ -189,7 +189,7 @@ end subroutine setComposition !! IT IS NOT THREAD SAFE! !! !! NOTE: eUpperSab and eLowerURR are fed by the aceNeutronDatabase, and they are the - !! strictest (respecitvely highest and lowest) energy limits among all nuclides + !! strictest (respectively highest and lowest) energy limits among all nuclides !! in the material composition. !! !! Args: @@ -219,7 +219,7 @@ subroutine set(self, name, matIdx, database, fissile, hasTMS, temp, eUpperSab, e if (present(fissile)) self % fissile = fissile if (present(matIdx)) self % matIdx = matIdx if (present(hasTMS)) self % hasTMS = hasTMS - if (present(temp)) self % kT = (kBoltzmann * temp) / joulesPerMeV + if (present(temp)) self % kT = kBoltzmannMeV * temp if (present(eUpperSab)) then if (eUpperSab < ZERO) call fatalError (Here, 'Upper Sab energy limit of material '& @@ -248,23 +248,30 @@ end subroutine set !! Args: !! xss [out] -> Cross section package to store the data !! E [in] -> Requested energy [MeV] + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> Random Number Generator !! !! Errors: !! fatalError if E is out-of-bounds for the stored data !! - subroutine getMacroXSs_byE(self, xss, E, rand) + subroutine getMacroXSs_byE(self, xss, E, temp, rho, rand) class(ceNeutronMaterial), intent(in) :: self type(neutronMacroXSs), intent(out) :: xss real(defReal), intent(in) :: E + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), intent(inout) :: rand - ! Check Cache and update if needed - if (materialCache(self % matIdx) % E_tail /= E .or. materialCache(self % matIdx) % E_tot /= E) then - call self % data % updateMacroXSs(E, self % matIdx, rand) - end if + associate( cache => materialCache(self % matIdx)) + ! Check Cache and update if needed + if (cache % E_tail /= E .or. cache % E_tot /= E .or.& + cache % T_tot /= temp .or. cache % rho_tot /= rho) then + call self % data % updateMacroXSs(E, self % matIdx, temp, rho, rand) + end if - xss = materialCache(self % matIdx) % xss + xss = cache % xss + end associate end subroutine getMacroXSs_byE @@ -321,51 +328,74 @@ end function useTMS !! rand [inout] -> random number generator !! nucIdx [out] -> sampled nuclide index !! eOut [out] -> relative energy between neutron and target (may be /= E in case of TMS) + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! !! Errors: !! fatalError if sampling fails for some reason (E.G. random number > 1) !! fatalError if E is out-of-bounds of the present data !! - subroutine sampleNuclide(self, E, rand, nucIdx, eOut) + subroutine sampleNuclide(self, E, rand, nucIdx, eOut, temp, rho) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E class(RNG), intent(inout) :: rand integer(shortInt), intent(out) :: nucIdx real(defReal), intent(out) :: eOut + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(ceNeutronNuclide), pointer :: nuc integer(shortInt) :: i real(defReal) :: P_acc, eMin, eMax, A, eRel, & - trackMatXS, totNucXS, dens + trackMatXS, totNucXS, dens, & + kT, densityFactor character(100), parameter :: Here = 'sampleNuclide (ceNeutronMaterial_class.f90)' ! Get material tracking XS - if (E /= materialCache(self % matIdx) % E_track) then - call self % data % updateTrackMatXS(E, self % matIdx, rand) + if (E /= materialCache(self % matIdx) % E_track .or. & + temp /= materialCache(self % matIdx) % T_track .or. & + rho /= materialCache(self % matIdx) % rho_track) then + call self % data % updateTrackMatXS(E, self % matIdx, temp, rho, rand) end if trackMatXS = materialCache(self % matIdx) % trackXS * rand % get() + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = self % kT + else + kT = temp * kBoltzmannMeV + end if + + ! Use imposed density if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if ! Loop over nuclides do i = 1,size(self % nuclides) nucIdx = self % nuclides(i) - dens = self % dens(i) + dens = self % dens(i) * densityFactor associate (nucCache => nuclideCache(nucIdx)) - + ! Retrieve nuclide XS from cache if (self % useTMS(E)) then ! If the material is using TMS, the nuclide temperature majorant is needed ! The check for the right values stored in cache happens inside the subroutine - call self % data % updateTotalTempNucXS(E, self % kT, nucIdx) + call self % data % updateTotalTempNucXS(E, kT, nucIdx) totNucXS = nucCache % tempMajXS * nucCache % doppCorr else ! Update nuclide cache if needed - if (E /= nucCache % E_tot) call self % data % updateTotalNucXS(E, nucIdx, self % kT, rand) + if (E /= nucCache % E_tot .or. temp /= materialCache(self % matIdx) % T_tot) then + call self % data % updateTotalNucXS(E, nucIdx, kT, rand) + end if totNucXS = nucCache % xss % total end if @@ -397,7 +427,7 @@ subroutine sampleNuclide(self, E, rand, nucIdx, eOut) if (eMax < eRel) eRel = eMax ! Get relative energy nuclide cross section - call self % data % updateTotalNucXS(eRel, nucIdx, self % kT, rand) + call self % data % updateTotalNucXS(eRel, nucIdx, kT, rand) ! Calculate acceptance probability using ratio of relative energy xs to temperature majorant P_acc = nucCache % xss % total * nucCache % doppCorr / totNucXS @@ -436,6 +466,8 @@ end subroutine sampleNuclide !! !! Args: !! E [in] -> incident energy [MeV] + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! !! Result: @@ -446,13 +478,15 @@ end subroutine sampleNuclide !! fatalError if E is out-of-bounds of the present data !! Returns nucIdx <= if material is not fissile !! - function sampleFission(self, E, rand) result(nucIdx) + function sampleFission(self, E, temp, rho, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), intent(inout) :: rand class(ceNeutronNuclide), pointer :: nuc integer(shortInt) :: nucIdx, i - real(defReal) :: xs, doppCorr, A, nuckT, deltakT + real(defReal) :: xs, doppCorr, A, nuckT, deltakT, kT, densityFactor character(100), parameter :: Here = 'sampleFission (ceNeutronMaterial_class.f90)' ! Short-cut for nonFissile material @@ -461,11 +495,25 @@ function sampleFission(self, E, rand) result(nucIdx) return end if + ! Used imposed density scaling if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = self % kT + else + kT = temp * kBoltzmannMeV + end if + ! Calculate material macroscopic nuFission ! The cache is updated without checking the energy to get the correct results with TMS ! The relative energy flag cached is cleaned to make sure cross sections are updated materialCache(self % matIdx) % E_rel = ZERO - call self % data % updateMacroXSs(E, self % matIdx, rand) + call self % data % updateMacroXSs(E, self % matIdx, temp, rho, rand) xs = materialCache(self % matIdx) % xss % nuFission * rand % get() @@ -482,7 +530,7 @@ function sampleFission(self, E, rand) result(nucIdx) A = nuc % getMass() nuckT = nuc % getkT() - deltakT = self % kT - nuckT + deltakT = kT - nuckT doppCorr = dopplerCorrectionFactor(E, A, deltakT) else @@ -492,7 +540,7 @@ function sampleFission(self, E, rand) result(nucIdx) ! The nuclide cache should be at the right energy after updating the material ! In the case of TMS where the macro xss are at a relative energy, the nuclide ! xss to be used are at the relative energy just sampled - xs = xs - nuclideCache(nucIdx) % xss % nuFission * self % dens(i) * doppCorr + xs = xs - nuclideCache(nucIdx) % xss % nuFission * self % dens(i) * doppCorr * densityFactor if (xs < ZERO) return @@ -511,6 +559,8 @@ end function sampleFission !! !! Args: !! E [in] -> incident energy [MeV] + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! !! Result: @@ -521,23 +571,39 @@ end function sampleFission !! fatalError if E is out-of-bounds of the present data !! Returns nucIdx <= if material is a pure-absorber (with fission as absorbtion) !! - function sampleScatter(self, E, rand) result(nucIdx) + function sampleScatter(self, E, temp, rho, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), intent(inout) :: rand class(ceNeutronNuclide), pointer :: nuc integer(shortInt) :: nucIdx, i - real(defReal) :: xs, doppCorr, A, nuckT, deltakT + real(defReal) :: xs, doppCorr, A, nuckT, deltakT, densityFactor, kT character(100), parameter :: Here = 'sampleScatter (ceNeutronMaterial_class.f90)' ! Calculate material macroscopic cross section of all scattering ! The cache is updated without checking the energy to get the correct results with TMS ! The relative energy flag cached is cleaned to make sure cross sections are updated materialCache(self % matIdx) % E_rel = ZERO - call self % data % updateMacroXSs(E, self % matIdx, rand) + call self % data % updateMacroXSs(E, self % matIdx, temp, rho, rand) xs = rand % get() * (materialCache(self % matIdx) % xss % elasticScatter + & materialCache(self % matIdx) % xss % inelasticScatter) + + ! Used imposed density scaling if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = self % kT + else + kT = temp * kBoltzmannMeV + end if ! Loop over all nuclides do i = 1,size(self % nuclides) @@ -553,7 +619,7 @@ function sampleScatter(self, E, rand) result(nucIdx) A = nuc % getMass() nuckT = nuc % getkT() - deltakT = self % kT - nuckT + deltakT = kT - nuckT doppCorr = dopplerCorrectionFactor(E, A, deltakT) else @@ -564,7 +630,8 @@ function sampleScatter(self, E, rand) result(nucIdx) ! In the case of TMS where the macro xss are at a relative energy, the nuclide ! xss to be used are at the relative energy just sampled xs = xs - (nuclideCache(nucIdx) % xss % elasticScatter + & - nuclideCache(nucIdx) % xss % inelasticScatter ) * self % dens(i) * doppCorr + nuclideCache(nucIdx) % xss % inelasticScatter ) & + * self % dens(i) * doppCorr * densityFactor if (xs < ZERO) return @@ -583,6 +650,8 @@ end function sampleScatter !! !! Args: !! E [in] -> incident energy [MeV] + !! temp [in] -> local temperature in kelvin. Negative values should be ignored. + !! rho [in] -> local density scaling factor. Negative values should be ignored. !! rand [inout] -> random number generator !! !! Result: @@ -593,25 +662,41 @@ end function sampleScatter !! fatalError if E is out-of-bounds of the present data !! Returns nucIdx <= if material is a pure-capture (with fission as scattering) !! - function sampleScatterWithFission(self, E, rand) result(nucIdx) + function sampleScatterWithFission(self, E, temp, rho, rand) result(nucIdx) class(ceNeutronMaterial), intent(in) :: self real(defReal), intent(in) :: E + real(defReal), intent(in) :: temp + real(defReal), intent(in) :: rho class(RNG), intent(inout) :: rand class(ceNeutronNuclide), pointer :: nuc integer(shortInt) :: nucIdx, i - real(defReal) :: xs, doppCorr, A, nuckT, deltakT + real(defReal) :: xs, doppCorr, A, nuckT, deltakT, kT, densityFactor character(100), parameter :: Here = 'sampleScatterWithFission (ceNeutronMaterial_class.f90)' ! Calculate material macroscopic cross section of all scattering and fission ! The cache is updated without checking the energy to get the correct results with TMS ! The relative energy flag cached is cleaned to make sure cross sections are updated materialCache(self % matIdx) % E_rel = ZERO - call self % data % updateMacroXSs(E, self % matIdx, rand) + call self % data % updateMacroXSs(E, self % matIdx, temp, rho, rand) xs = rand % get() * (materialCache(self % matIdx) % xss % elasticScatter + & materialCache(self % matIdx) % xss % inelasticScatter + & materialCache(self % matIdx) % xss % fission) + ! Used imposed density scaling if given + if (rho <= ZERO) then + densityFactor = ONE + else + densityFactor = rho + end if + + ! Use imposed temperature if given + if (temp <= ZERO) then + kT = self % kT + else + kT = temp * kBoltzmannMeV + end if + ! Loop over all nuclides do i = 1,size(self % nuclides) @@ -626,7 +711,7 @@ function sampleScatterWithFission(self, E, rand) result(nucIdx) A = nuc % getMass() nuckT = nuc % getkT() - deltakT = self % kT - nuckT + deltakT = kT - nuckT doppCorr = dopplerCorrectionFactor(E, A, deltakT) else @@ -638,7 +723,8 @@ function sampleScatterWithFission(self, E, rand) result(nucIdx) ! xss to be used are at the relative energy just sampled xs = xs - (nuclideCache(nucIdx) % xss % elasticScatter + & nuclideCache(nucIdx) % xss % inelasticScatter + & - nuclideCache(nucIdx) % xss % fission) * self % dens(i) * doppCorr + nuclideCache(nucIdx) % xss % fission) * & + self % dens(i) * doppCorr * densityFactor if (xs < ZERO) return diff --git a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 index cf824d6d3..d874eca34 100644 --- a/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 +++ b/NuclearData/mgNeutronData/baseMgNeutron/baseMgNeutronDatabase_class.f90 @@ -74,9 +74,9 @@ module baseMgNeutronDatabase_class procedure :: kill procedure :: init procedure :: activate + procedure :: initMajorant ! Local interface - procedure :: initMajorant procedure :: nGroups end type baseMgNeutronDatabase @@ -407,7 +407,7 @@ subroutine init(self, dict, ptr, silent) self % nG = self % mats(1) % nGroups() do i = 2,nMat if(self % nG /= self % mats(i) % nGroups()) then - call fatalError(Here,'Inconsistant # of groups in materials in matIdx'//numToChar(i)) + call fatalError(Here,'Inconsistent # of groups in materials in matIdx '//numToChar(i)) end if end do @@ -455,13 +455,36 @@ end subroutine activate !! !! See nuclearDatabase documentation for details !! - subroutine initMajorant(self, loud) + subroutine initMajorant(self, loud, maxTemp, scaleDensity) class(baseMgNeutronDatabase), intent(inout) :: self - logical(defBool), intent(in) :: loud + logical(defBool), optional, intent(in) :: loud + real(defReal), optional, intent(in) :: maxTemp + real(defReal), optional, intent(in) :: scaleDensity + logical(defBool) :: isLoud integer(shortInt) :: g, i, idx - real(defReal) :: xs + real(defReal) :: xs, densityFactor integer(shortInt), parameter :: TOTAL_XS = 1 + if (present(loud)) then + isLoud = loud + else + isLoud = .false. + end if + + ! Scale density + if (present(scaleDensity)) then + if (scaleDensity < ONE) then + densityFactor = ONE + else + densityFactor = scaleDensity + end if + else + densityFactor = ONE + end if + + ! Currently ignores maxTemp input + ! TODO: Update should there be a temperature model developed for MG XSs + ! Allocate majorant allocate (self % majorant(self % nG)) @@ -472,10 +495,10 @@ subroutine initMajorant(self, loud) idx = self % activeMats(i) xs = max(xs, self % mats(idx) % data(TOTAL_XS, g)) end do - self % majorant(g) = xs + self % majorant(g) = xs * densityFactor end do - if (loud) print '(A)', 'MG unionised majorant cross section calculation completed' + if (isLoud) print '(A)', 'MG unionised majorant cross section calculation completed' end subroutine initMajorant diff --git a/NuclearData/nuclearDatabase_inter.f90 b/NuclearData/nuclearDatabase_inter.f90 index bada8a797..54acc80c0 100644 --- a/NuclearData/nuclearDatabase_inter.f90 +++ b/NuclearData/nuclearDatabase_inter.f90 @@ -40,6 +40,7 @@ module nuclearDatabase_inter contains procedure(init), deferred :: init procedure(activate), deferred :: activate + procedure(initMajorant), deferred :: initMajorant procedure(getTrackingXS), deferred :: getTrackingXS procedure(getTrackMatXS), deferred :: getTrackMatXS procedure(getTotalMatXS), deferred :: getTotalMatXS @@ -90,12 +91,39 @@ subroutine activate(self, activeMat, silent) logical(defBool), optional, intent(in) :: silent end subroutine activate + !! + !! Constructs the majorant cross section. + !! Can be called repeatedly to update the majorant subject + !! to changes in the geometry. + !! + !! Optionally can be 'loud', i.e., outputs information on initialising + !! the majorant. + !! + !! Optionally can account for the temperature due to super-imposed + !! temperature fields. This should receive the maximum temperature in + !! the system. This is conservative and can be improved by creating a + !! material-wise maximum temperature input. + !! + !! Optionally can scale the density. scaleDensity should be + !! the relative (to input) density of the highest density material. + !! This is most naturally used with a super-imposed density field. + !! As for temperature, this approach for density is conservative and + !! could be improved with a material-wise density scaling factor. + !! + subroutine initMajorant(self, loud, maxTemp, scaleDensity) + import :: nuclearDatabase, defBool, defReal + class(nuclearDatabase), intent(inout) :: self + logical(defBool), optional, intent(in) :: loud + real(defReal), optional, intent(in) :: maxTemp + real(defReal), optional, intent(in) :: scaleDensity + end subroutine initMajorant + !! !! Return value of tracking XS for a particle and a given request !! !! Reads all relevant state information from the particle (e.g. E or G) !! It is the XS used to sample track length: it might be the same as the - !! material total XS, the majorant XS, or a temperature majorant + !! material total XS, the majorant XS, or a temperature majorant. !! !! Args: !! p [in] -> Particle at a given state @@ -106,8 +134,9 @@ end subroutine activate !! Value of XS used to sample path length [1/cm] !! !! Errors: - !! Undefined behaviour if the state of the particle is invalid e.g. -ve energy - !! Undefined behavior if matIdx does not correspond to a defined material + !! Undefined behaviour if the state of the particle is invalid e.g. -ve energy. + !! Undefined behavior if matIdx does not correspond to a defined material. + !! Fatal error if incorrect 'what' option is provided. !! function getTrackingXS(self, p, matIdx, what) result(xs) import :: nuclearDatabase, particle, shortInt, defReal diff --git a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 index e6032e73b..32f17eb62 100644 --- a/NuclearData/testNeutronData/testNeutronDatabase_class.f90 +++ b/NuclearData/testNeutronData/testNeutronDatabase_class.f90 @@ -44,6 +44,7 @@ module testNeutronDatabase_class ! Superclass Interface procedure :: init procedure :: activate + procedure :: initMajorant procedure :: getTrackingXS procedure :: getTrackMatXS procedure :: getTotalMatXS @@ -162,6 +163,19 @@ subroutine activate(self, activeMat, silent) ! Do nothing end subroutine activate + + !! + !! + !! + subroutine initMajorant(self, loud, maxTemp, scaleDensity) + class(testNeutronDatabase), intent(inout) :: self + logical(defBool), intent(in), optional :: loud + real(defReal), intent(in), optional :: maxTemp + real(defReal), intent(in), optional :: scaleDensity + + ! Do nothing + + end subroutine initMajorant !! !! Return value of Tracking XS for a particle and a given request diff --git a/ParticleObjects/Source/fissionSource_class.f90 b/ParticleObjects/Source/fissionSource_class.f90 index 3c2da3ca3..b9117af46 100644 --- a/ParticleObjects/Source/fissionSource_class.f90 +++ b/ParticleObjects/Source/fissionSource_class.f90 @@ -2,7 +2,8 @@ module fissionSource_class use numPrecision use endfConstants - use universalVariables, only : OUTSIDE_MAT, VOID_MAT, UNDEF_MAT, OVERLAP_MAT + use universalVariables, only : OUTSIDE_MAT, VOID_MAT, UNDEF_MAT, OVERLAP_MAT, & + NO_TEMPERATURE, NO_DENSITY use genericProcedures, only : rotateVector, numToChar use errors_mod, only : fatalError use dictionary_class, only : dictionary @@ -224,7 +225,7 @@ function sampleParticle(self, rand) result(p) end if ! Get Nuclide - nucIdx = matCE % sampleFission(self % E, rand) + nucIdx = matCE % sampleFission(self % E, NO_TEMPERATURE, NO_DENSITY, rand) ! Get reaction object fissCE => fissionCE_TptrCast(nucData % getReaction(N_FISSION, nucIdx)) diff --git a/ParticleObjects/Source/materialSource_class.f90 b/ParticleObjects/Source/materialSource_class.f90 index 37438c806..173b0e549 100644 --- a/ParticleObjects/Source/materialSource_class.f90 +++ b/ParticleObjects/Source/materialSource_class.f90 @@ -187,7 +187,7 @@ function sampleParticle(self, rand) result(p) call self % geom % whatIsAt(matIdx, uniqueID, r) ! Reject if there is no material or if the particle is in void - if (matIdx == OUTSIDE_MAT .or. matIdx == VOID_MAT) cycle rejection + if (matIdx == OUTSIDE_MAT) cycle rejection mat => neutronMaterial_CptrCast(nucData % getMaterial(matIdx)) if (.not.associated(mat)) call fatalError(Here, "Nuclear data did not return neutron material.") diff --git a/ParticleObjects/particle_class.f90 b/ParticleObjects/particle_class.f90 index 0e06a4115..818033f11 100644 --- a/ParticleObjects/particle_class.f90 +++ b/ParticleObjects/particle_class.f90 @@ -111,6 +111,10 @@ module particle_class real(defReal) :: w ! Particle Weight real(defReal) :: time ! Particle time point + ! Information passed from geometry + real(defReal) :: T = NO_TEMPERATURE ! Local temperature + real(defReal) :: rho = NO_DENSITY ! Local density scaling + ! Precursor particle data real(defReal) :: lambda = INF ! Precursor decay constant diff --git a/PhysicsPackages/eigenPhysicsPackage_class.f90 b/PhysicsPackages/eigenPhysicsPackage_class.f90 index 878780e4f..76829770e 100644 --- a/PhysicsPackages/eigenPhysicsPackage_class.f90 +++ b/PhysicsPackages/eigenPhysicsPackage_class.f90 @@ -31,11 +31,13 @@ module eigenPhysicsPackage_class ! Geometry use geometry_inter, only : geometry use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx, & - gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr + gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_fieldPtrName => fieldPtrName use geometryFactory_func, only : new_geometry ! Fields use field_inter, only : field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast use uniFissSitesField_class, only : uniFissSitesField, uniFissSitesField_TptrCast use fieldFactory_func, only : new_field @@ -426,6 +428,8 @@ subroutine init(self, dict) type(outputFile) :: test_out type(visualiser) :: viz class(field), pointer :: field + class(pieceConstantField), pointer :: pcField + real(defReal) :: maxDensityScale, maxTemperature character(100), parameter :: Here ='init (eigenPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -522,6 +526,31 @@ subroutine init(self, dict) call viz % makeViz() call viz % kill() endif + + ! If present, build temperature field + if (dict % isPresent('temperature')) then + tempDict => dict % getDictPtr('temperature') + call new_field(tempDict, nameTemperature) + field => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(field) + maxTemperature = pcField % getMaxValue() + else + maxTemperature = NO_TEMPERATURE + end if + + ! If present, build density field + if (dict % isPresent('density')) then + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + field => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(field) + maxDensityScale = pcField % getMaxValue() + else + maxDensityScale = NO_DENSITY + end if + + ! Update majorant in case of density and temperature fields + call self % nucData % initMajorant(.false., maxTemp = maxTemperature, scaleDensity = maxDensityScale) ! Read uniform fission site option as a geometry field if (dict % isPresent('uniformFissionSites')) then diff --git a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 index e5e664591..1489b0f08 100644 --- a/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 +++ b/PhysicsPackages/fixedSourcePhysicsPackage_class.f90 @@ -30,11 +30,15 @@ module fixedSourcePhysicsPackage_class ! Geometry use geometry_inter, only : geometry - use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx + use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx, & + gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_fieldPtrName => fieldPtrName use geometryFactory_func, only : new_geometry ! Fields use fieldFactory_func, only : new_field + use field_inter, only : field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast ! Nuclear Data use materialMenu_mod, only : mm_nMat => nMat @@ -331,6 +335,9 @@ subroutine init(self, dict) character(nameLen) :: nucData, energy, geomName type(outputFile) :: test_out type(visualiser) :: viz + class(field), pointer :: field + class(pieceConstantField), pointer :: pcField + real(defReal) :: maxTemperature, maxDensityScale character(100), parameter :: Here ='init (fixedSourcePhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -408,6 +415,31 @@ subroutine init(self, dict) ! Activate Nuclear Data *** All materials are active call ndReg_activate(self % particleType, nucData, self % geom % activeMats()) self % nucData => ndReg_get(self % particleType) + + ! If present, build temperature field + if (dict % isPresent('temperature')) then + tempDict => dict % getDictPtr('temperature') + call new_field(tempDict, nameTemperature) + field => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(field) + maxTemperature = pcField % getMaxValue() + else + maxTemperature = NO_TEMPERATURE + end if + + ! If present, build density field + if (dict % isPresent('density')) then + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + field => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(field) + maxDensityScale = pcField % getMaxValue() + else + maxDensityScale = NO_DENSITY + end if + + ! Update majorant in case of density and temperature fields + call self % nucData % initMajorant(.false., maxTemp = maxTemperature, scaleDensity = maxDensityScale) ! Call visualisation if (dict % isPresent('viz') .and. isMPIMaster()) then diff --git a/PhysicsPackages/kineticPhysicsPackage_class.f90 b/PhysicsPackages/kineticPhysicsPackage_class.f90 index 9ac88459d..4ca5eeac5 100644 --- a/PhysicsPackages/kineticPhysicsPackage_class.f90 +++ b/PhysicsPackages/kineticPhysicsPackage_class.f90 @@ -29,8 +29,15 @@ module kineticPhysicsPackage_class ! Geometry use geometry_inter, only : geometry - use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx + use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx, & + gr_fieldIdx => fieldIdx, gr_fieldPtr => fieldPtr, & + gr_fieldPtrName => fieldPtrName use geometryFactory_func, only : new_geometry + + ! Fields + use field_inter, only : field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast + use fieldFactory_func, only : new_field ! Nuclear Data use materialMenu_mod, only : mm_nMat => nMat @@ -498,6 +505,9 @@ subroutine init(self, dict) character(nameLen) :: nucData, energy, geomName type(outputFile) :: test_out type(visualiser) :: viz + class(field), pointer :: field + class(pieceConstantField), pointer :: pcField + real(defReal) :: maxTemperature, maxDensityScale character(100), parameter :: Here ='init (kineticPhysicsPackage_class.f90)' call cpu_time(self % CPU_time_start) @@ -618,6 +628,31 @@ subroutine init(self, dict) ! Activate Nuclear Data *** All materials are active call ndReg_activate(self % particleType, nucData, self % geom % activeMats()) self % nucData => ndReg_get(self % particleType) + + ! If present, build temperature field + if (dict % isPresent('temperature')) then + tempDict => dict % getDictPtr('temperature') + call new_field(tempDict, nameTemperature) + field => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(field) + maxTemperature = pcField % getMaxValue() + else + maxTemperature = NO_TEMPERATURE + end if + + ! If present, build density field + if (dict % isPresent('density')) then + tempDict => dict % getDictPtr('density') + call new_field(tempDict, nameDensity) + field => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(field) + maxDensityScale = pcField % getMaxValue() + else + maxDensityScale = NO_DENSITY + end if + + ! Update majorant in case of density and temperature fields + call self % nucData % initMajorant(.false., maxTemp = maxTemperature, scaleDensity = maxDensityScale) ! Call visualisation if (dict % isPresent('viz') .and. isMPIMaster()) then diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 790fff3e5..7986d4fc3 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -28,7 +28,8 @@ module universalVariables integer(shortInt), parameter, public :: COLL_EV = 1, & BOUNDARY_EV = 2, & CROSS_EV = 3, & - LOST_EV = 4 + LOST_EV = 4, & + FIELD_EV = 5 ! Create definitions for readability when dealing with positions relative to surfaces logical(defBool), parameter, public :: behind = .FALSE., & @@ -75,22 +76,29 @@ module universalVariables integer(shortInt), parameter :: MATERIAL_XS = 1, & MAJORANT_XS = 2, & TRACKING_XS = 3 + + ! Unit conversion + real(defReal), parameter :: joulesPerMeV = 1.60218e-13_defReal ,& ! Convert MeV to J + shakesPerS = 1.0e+8_defReal ! Convert shakes to s ! Physical constants ! Neutron mass and speed of light in vacuum from from https://physics.nist.gov/cuu/Constants/index.html real(defReal), parameter :: neutronMass = 939.56542194_defReal, & ! Neutron mass in MeV (m*c^2) lightSpeed = 2.99792458e10_defReal, & ! Light speed in cm/s kBoltzmann = 1.380649e-23_defReal, & ! Bolztmann constant in J/K - energyPerFission = 200.0_defReal ! MeV - - ! Unit conversion - real(defReal), parameter :: joulesPerMeV = 1.60218e-13_defReal ,& ! Convert MeV to J - shakesPerS = 1.0e8_defReal ! Convert shakes to s + energyPerFission = 200.0_defReal, & ! MeV + kBoltzmannMeV = kBoltzmann / joulesPerMeV ! Global name variables used to define specific geometry or field types - character(nameLen), parameter :: nameUFS = 'uniFissSites' - character(nameLen), parameter :: nameWW = 'WeightWindows' - + character(nameLen), parameter :: nameUFS = 'uniFissSites' + character(nameLen), parameter :: nameWW = 'WeightWindows' + character(nameLen), parameter :: nameTemperature = 'temperature' + character(nameLen), parameter :: nameDensity = 'density' + + ! Flags associated with fields + real(defReal), parameter :: NO_TEMPERATURE = -INF, & + NO_DENSITY = -INF + ! Flag to indicate source file format integer(shortInt), parameter, public :: NO_PRINTING = 0, & ASCII_FILE = 1, & diff --git a/TransportOperator/transportOperatorDT_class.f90 b/TransportOperator/transportOperatorDT_class.f90 index 6d0d6d067..8aefa39d7 100644 --- a/TransportOperator/transportOperatorDT_class.f90 +++ b/TransportOperator/transportOperatorDT_class.f90 @@ -5,25 +5,25 @@ module transportOperatorDT_class use numPrecision use universalVariables - use errors_mod, only : fatalError - use genericProcedures, only : numToChar - use particle_class, only : particle - use particleDungeon_class, only : particleDungeon - use dictionary_class, only : dictionary + use errors_mod, only : fatalError + use genericProcedures, only : numToChar + use particle_class, only : particle + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary ! Superclass - use transportOperator_inter, only : transportOperator, init_super => init + use transportOperator_inter, only : transportOperator, init_super => init ! Geometry interfaces - use geometry_inter, only : geometry + use geometry_inter, only : geometry ! Tally interface use tallyCodes - use tallyAdmin_class, only : tallyAdmin + use tallyAdmin_class, only : tallyAdmin ! Nuclear data interfaces - use nuclearDataReg_mod, only : ndReg_get => get - use nuclearDatabase_inter, only : nuclearDatabase + use nuclearDataReg_mod, only : ndReg_get => get + use nuclearDatabase_inter, only : nuclearDatabase implicit none private @@ -109,6 +109,9 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) exit DTLoop end if + ! Get local conditions of temperature and density + call self % localConditions(p) + ! Obtain the local cross-section sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index 2c505f489..c1461b4a0 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -62,6 +62,9 @@ subroutine tracking_selection(self, p, tally, thisCycle, nextCycle) if (p % matIdx() == VOID_MAT) then sigmaT = ZERO else + ! Get local conditions + call self % localConditions(p) + sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) end if @@ -145,6 +148,9 @@ subroutine deltaTracking(self, p, tally, thisCycle, nextCycle) exit DTLoop end if + ! Get local conditions + call self % localConditions(p) + ! Obtain the local cross-section sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) @@ -180,6 +186,9 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) STLoop: do + ! Get local conditions + call self % localConditions(p) + sigmaTrack = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS) ! Obtain the local cross-section, depending on the material diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index 3b287b3a3..f73d27563 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -5,23 +5,23 @@ module transportOperatorST_class use numPrecision use universalVariables - use errors_mod, only : fatalError - use particle_class, only : particle - use particleDungeon_class, only : particleDungeon - use dictionary_class, only : dictionary + use errors_mod, only : fatalError + use particle_class, only : particle + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary ! Superclass - use transportOperator_inter, only : transportOperator, init_super => init + use transportOperator_inter, only : transportOperator, init_super => init ! Geometry interfaces - use geometry_inter, only : geometry, distCache + use geometry_inter, only : geometry, distCache ! Tally interface use tallyCodes - use tallyAdmin_class, only : tallyAdmin + use tallyAdmin_class, only : tallyAdmin ! Nuclear data interfaces - use nuclearDatabase_inter, only : nuclearDatabase + use nuclearDatabase_inter, only : nuclearDatabase implicit none private @@ -61,6 +61,9 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) STLoop: do + ! Get local conditions + call self % localConditions(p) + sigmaTrack = self % xsData % getTrackingXS(p, p % matIdx(), MATERIAL_XS) ! Obtain the local cross-section, depending on the material @@ -72,7 +75,7 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) sigmaT = ZERO else - + invSigmaTrack = ONE / sigmaTrack dist = -log( p % pRNG % get()) * invSigmaTrack @@ -80,7 +83,12 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) sigmaT = self % xsData % getTrackMatXS(p, p % matIdx()) ! Should never happen! Catches NaN distances - if (dist /= dist) call fatalError(Here, "Distance is NaN") + if (dist /= dist) then + print *, "Particle location: ", p % rGlobal() + print *, "Particle direction: ", p % dirGlobal() + print *, "Total XS: ", sigmaT + call fatalError(Here, "Distance is NaN") + end if end if diff --git a/TransportOperator/transportOperator_inter.f90 b/TransportOperator/transportOperator_inter.f90 index 6e035b24a..3659a48f6 100644 --- a/TransportOperator/transportOperator_inter.f90 +++ b/TransportOperator/transportOperator_inter.f90 @@ -2,23 +2,28 @@ module transportOperator_inter use numPrecision use universalVariables - use errors_mod, only : fatalError + use errors_mod, only : fatalError - use particle_class, only : particle - use particleDungeon_class, only : particleDungeon - use dictionary_class, only : dictionary + use particle_class, only : particle + use particleDungeon_class, only : particleDungeon + use dictionary_class, only : dictionary ! Geometry interfaces - use geometryReg_mod, only : gr_geomPtr => geomPtr - use geometry_inter, only : geometry + use geometryReg_mod, only : gr_geomPtr => geomPtr, & + gr_hasField => hasField, & + gr_fieldPtrName => fieldPtrName + use geometry_inter, only : geometry + + use field_inter, only : field + use pieceConstantField_inter, only : pieceConstantField, pieceConstantField_CptrCast ! Tally interface - use tallyAdmin_class, only : tallyAdmin + use tallyAdmin_class, only : tallyAdmin ! Nuclear data interfaces - use nuclearDataReg_mod, only : ndReg_get => get - use nuclearDatabase_inter, only : nuclearDatabase + use nuclearDataReg_mod, only : ndReg_get => get + use nuclearDatabase_inter, only : nuclearDatabase @@ -42,6 +47,9 @@ module transportOperator_inter !! Customisable procedures or transport actions !! transit(p, tally, thisCycle, nextCycle) -> implements movement from collision to collision !! + !! Procedures generic to transport operators: + !! localConditions(p) -> obtains local conditions of temperature and density for use in transport + !! type, abstract, public :: transportOperator !! Nuclear Data block pointer -> public so it can be used by subclasses (protected member) class(nuclearDatabase), pointer :: xsData => null() @@ -57,6 +65,9 @@ module transportOperator_inter procedure :: init procedure :: kill + ! Query for local conditions of temperature and density. + procedure, non_overridable :: localConditions + ! Customisable deferred procedures procedure(transit), deferred :: transit @@ -118,6 +129,32 @@ subroutine transport(self, p, tally, thisCycle, nextCycle) end subroutine transport + !! + !! Queries local conditions of temperature and density from + !! the geometry registry. Updates the particle to carry this info. + !! + subroutine localConditions(self, p) + class(transportOperator), intent(in) :: self + class(particle), intent(inout) :: p + class(field), pointer :: genericField + class(pieceConstantField), pointer :: pcField + + ! Temperature check + if (gr_hasField(nameTemperature)) then + genericField => gr_fieldPtrName(nameTemperature) + pcField => pieceConstantField_CptrCast(genericField) + p % T = pcField % at(p % coords) + end if + + ! Density check + if (gr_hasField(nameDensity)) then + genericField => gr_fieldPtrName(nameDensity) + pcField => pieceConstantField_CptrCast(genericField) + p % rho = pcField % at(p % coords) + end if + + end subroutine localConditions + !! !! Initialise transport operator from dictionary and geometry !! diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index 9416b1a23..f3d214b40 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -593,6 +593,8 @@ A detailed description about the geometry modelling adopted in SCONE can be foun surfaces { } cells { } universes { } + temperature { } + density { } } At the moment, the only **geometry** type available is ``geometryStd``. As for the boundary @@ -628,6 +630,12 @@ Hence, an example of a geometry input could look like: :: For more details about the graph-like structure of the nested geometry see the relevant :ref:`section `. +The geometry optionally allows the use of ``temperature`` and ``density`` fields. These +are super-imposed fields which modify the temperature and densities given in nuclear data. +The temperature field specifies local temperatures in kelvin while the density field +specifies the local dimensionless factors by which material density should be scaled. +Both of these follow the syntax of a ``PieceConstantField``. + Surfaces ######## @@ -902,7 +910,7 @@ Example: :: 1 2 3 // x: 1-3, y: 2, z: 2 4 5 6 // x: 1-3, y: 1, z: 2 7 8 9 // x: 1-3, y: 2, z: 1 - 10 11 12 ) } // x: 1-3, y: 1, z: 1 + 10 11 12 ); } // x: 1-3, y: 1, z: 1 .. note:: The order of the elements in the lattice is different from other MC codes, e.g., @@ -917,6 +925,52 @@ Example: :: root { id 1000; type rootUniverse; border 10; fill u<1>; } +PieceConstantFields +################### + +These are fields which are piecewise constant and are endowed with a distance calculation to +compute the distance until the value of the field changes. These can be used for imposing +density and temperature distributions across the system in a convenient manner. Can be initialised +either with an explicit definition or with a path to the field definition. + +Currently there is only one available PieceConstantField: + +* cartesianField. This is similar to a latUniverse: the value of the field varies over a regular + Cartesian lattice with a given shape and size. The field also allows specifying different values + in different materials, or uniformly across all materials. + + - shape: (x y z) array of integers, stating the numbers of x, y and z + elements of the field. For a 2D field, the z entry has to be 0. + - pitch: (x y z) array with the x, y and z field pitches. In a 2D field, + the value entered in the third dimension is not used. [cm] + - origin (*optional*, default = (0.0 0.0 0.0)): (x y z) array with the + origin of the field. [cm] + - materials: list of material names, corresponding to materials in nuclearData. + Optionally, ``all`` can be used, applying the values of the field to all materials. + - names of each material: a map, named after every material present in the materials list. + The entries of the map are the values that the field takes in that material in that + element of the field. The order is: increasing x, increasing y and then increasing z. + - default: the value taken by the field when a point is either outside of the field or + in a material which is not included in the field. + +Example: :: + + temperature { type cartesianField; shape (3 2 2); pitch (1.0 1.0 1.5); + materials (uo2 water); + uo2 ( + 901 902 903 + 904 905 906 + 907 908 909 + 910 911 912 ); + water ( + 601 602 603 + 604 605 606 + 607 608 609 + 610 611 612); + default 302; } + + density { type cartesianField; file ./myDensityField; } + Visualiser ---------- @@ -1634,6 +1688,14 @@ Examples: :: map1 { type spaceMap; axis x; grid lin; min -50.0; max 50.0; N 100; } map2 { type spaceMap; axis z; grid unstruct; bins (0.0 0.2 0.3 0.5 0.7 0.8 1.0); } +* fieldMap (1D map), map over superimposed fields. Limited currently to pieceConstantFields. + + - field: field definition, corresponding to those in pieceConstantFields. + +Examples: :: + + map1 { type fieldMap; field {file ./myField.txt } } + * timeMap (1D map), maps particles point in time - grid: ``lin`` for linearly spaced bins