Temperature and density fields

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

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()

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()

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,22 +347,21 @@ elemental subroutine kill(self)
     self % nMat = 0
     self % corner = ZERO
     self % a_bar  = ZERO
-    self % outLocalID = 0
     call self % outline % kill()
 
   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))