Implicit MC

CMakeLists.txt

@@ -77,6 +77,9 @@ list(APPEND CMAKE_PREFIX_PATH $ENV{LAPACK_INSTALL})
 find_package(LAPACK REQUIRED )
 message(STATUS ${LAPACK_LIBRARIES})
 
+set(BLAS_LIBS /home/ajb343/BLAS/BLAS-3.10.0/Build)
+find_package(BLAS REQUIRED)
+
 # Dependencies for BUILD_TESTS
 if (BUILD_TESTS)
     # FIND PYTHON INTERPRETER

CollisionOperator/CollisionProcessors/CMakeLists.txt

@@ -2,5 +2,6 @@
 add_sources( ./collisionProcessor_inter.f90 
              ./collisionProcessorFactory_func.f90 
              ./neutronCEstd_class.f90
-	     ./neutronCEimp_class.f90
-             ./neutronMGstd_class.f90) 
+             ./neutronCEimp_class.f90
+             ./neutronMGstd_class.f90
+             ./IMCMGstd_class.f90) 

CollisionOperator/CollisionProcessors/IMCMGstd_class.f90

@@ -0,0 +1,220 @@
+module IMCMGstd_class
+
+  use numPrecision
+  use endfConstants
+  use genericProcedures,             only : fatalError, rotateVector, numToChar
+  use dictionary_class,              only : dictionary
+  use RNG_class,                     only : RNG
+
+  ! Particle types
+  use particle_class,                only : particle, particleState, printType, P_PHOTON
+  use particleDungeon_class,         only : particleDungeon
+
+  ! Abstract interface
+  use collisionProcessor_inter,      only : collisionProcessor, collisionData ,init_super => init
+
+  ! Nuclear Data Interface
+  use nuclearDataReg_mod,            only : ndReg_getIMCMG => getIMCMG
+  use nuclearDatabase_inter,         only : nuclearDatabase
+  use mgIMCDatabase_inter,           only : mgIMCDatabase
+  use mgIMCMaterial_inter,           only : mgIMCMaterial, mgIMCMaterial_CptrCast
+  use reactionHandle_inter,          only : reactionHandle
+
+  implicit none
+  private
+
+  !!
+  !! Standard (default) scalar collision processor for MG IMC
+  !! Determines type of collision as either absorption or effective scattering
+  !!
+  !! Settings:
+  !!  NONE
+  !!
+  !! Sample dictionary input:
+  !!   collProcName {
+  !!   type            IMCMGstd;
+  !!   }
+  !!
+  type, public, extends(collisionProcessor) :: IMCMGstd
+    private
+    class(mgIMCDatabase), pointer, public :: xsData => null()
+    class(mgIMCMaterial), pointer, public :: mat    => null()
+  contains
+    ! Initialisation procedure
+    procedure :: init
+
+    ! Implementation of customisable procedures
+    procedure :: sampleCollision
+    procedure :: implicit
+    procedure :: elastic
+    procedure :: inelastic
+    procedure :: capture
+    procedure :: fission
+    procedure :: cutoffs
+  end type IMCMGstd
+
+contains
+
+  !!
+  !! Initialise from dictionary
+  !!
+  subroutine init(self, dict)
+    class(IMCMGstd), intent(inout)     :: self
+    class(dictionary), intent(in)      :: dict
+    character(100), parameter :: Here = 'init (IMCMGstd_class.f90)'
+
+    ! Call superclass
+    call init_super(self, dict)
+
+  end subroutine init
+
+  !!
+  !! Samples collision
+  !!
+  !! Absorption with probability equal to fleck factor, otherwise
+  !!  effective scattering
+  !!
+  !! Physical scattering is omitted as in reference paper "Four Decades of Implicit Monte Carlo"
+  !!  (Allan B Wollaber) but may be included later if desired 
+  !!
+  subroutine sampleCollision(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+    real(defReal)                        :: r, fleck
+    character(100),parameter :: Here =' sampleCollision (IMCMGstd_class.f90)'
+
+    ! Verify that particle is MG PHOTON
+    if( .not. p % isMG .or. p % type /= P_PHOTON) then
+      call fatalError(Here, 'Supports only MG PHOTON. Was given NEUTRON and/or CE '//printType(p % type))
+    end if
+
+    ! Verify and load nuclear data pointer
+    self % xsData => ndReg_getIMCMG()
+    if(.not.associated(self % xsData)) call fatalError(Here, "Failed to get active database for MG IMC")
+
+    ! Get and verify material pointer
+    self % mat => mgIMCMaterial_CptrCast( self % xsData % getMaterial( p % matIdx()))
+    if(.not.associated(self % mat)) call fatalError(Here, "Failed to get MG IMC Material")
+
+    r = p % pRNG % get()
+
+    fleck = self % mat % getFleck()
+
+    if( r < fleck ) then
+      ! Effective absoprtion
+      collDat % MT = macroCapture
+    else
+      ! Effective scattering
+      collDat % MT = macroAllScatter
+    end if
+
+  end subroutine sampleCollision
+
+  !!
+  !! Perform implicit treatment
+  !!
+  subroutine implicit(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+
+    ! Do nothing.
+
+  end subroutine implicit
+
+  !!
+  !! Effective scattering - currently only elastic (constant energy-weight)
+  !!
+  subroutine elastic(self, p , collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+    real(defReal)                        :: phi, mu
+    real(defReal), dimension(3)          :: dir
+    character(100), parameter :: Here = 'elastic (IMCMGstd_class.f90)'
+
+    ! Assign MT number
+    collDat % MT = macroAllScatter 
+
+    ! Sample Direction - chosen uniformly inside unit sphere
+    mu = 2 * p % pRNG % get() - 1
+    phi = p % pRNG % get() * 2*pi
+    dir(1) = mu
+    dir(2) = sqrt(1-mu**2) * cos(phi)
+    dir(3) = sqrt(1-mu**2) * sin(phi)
+
+    !p % coords % dir = dir
+    call p % rotate(mu, phi)
+
+  end subroutine elastic
+
+  !!
+  !! Inelastic scattering - Not currently supported
+  !!
+  subroutine inelastic(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+    character(100),parameter :: Here = "inelastic (IMCMGstd_class.f90)"
+
+    ! Do nothing. Should not be called
+
+    call fatalError(Here, "Inelastic subroutine should not be called")
+
+  end subroutine inelastic
+
+  !!
+  !! Perform capture
+  !!
+  subroutine capture(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+
+    p % isDead = .true.
+
+  end subroutine capture
+
+  !!
+  !! Perform fission
+  !!
+  subroutine fission(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+    character(100), parameter :: Here = 'fission (IMCMGstd_class.f90)'
+
+    ! Do nothing. Should not be called
+
+    call fatalError(Here, "Fission subroutine should not be called")
+
+  end subroutine fission
+
+  !!
+  !! Apply cutoffs or post-collision implicit treatment
+  !!
+  subroutine cutoffs(self, p, collDat, thisCycle, nextCycle)
+    class(IMCMGstd), intent(inout)       :: self
+    class(particle), intent(inout)       :: p
+    type(collisionData), intent(inout)   :: collDat
+    class(particleDungeon),intent(inout) :: thisCycle
+    class(particleDungeon),intent(inout) :: nextCycle
+
+    ! Do nothing
+
+  end subroutine cutoffs
+
+end module IMCMGstd_class

CollisionOperator/CollisionProcessors/collisionProcessorFactory_func.f90

@@ -11,6 +11,7 @@ module collisionProcessorFactory_func
   use neutronCEstd_class, only : neutronCEstd
   use neutronCEimp_class, only : neutronCEimp
   use neutronMGstd_class, only : neutronMGstd
+  use IMCMGstd_class,     only : IMCMGstd
 
   implicit none
   private
@@ -24,7 +25,8 @@ module collisionProcessorFactory_func
   ! For now  it is necessary to adjust trailing blanks so all enteries have the same length
   character(nameLen),dimension(*),parameter :: AVALIBLE_collisionProcessors = [ 'neutronCEstd',&
                                                                                 'neutronCEimp',&
-                                                                                'neutronMGstd']
+                                                                                'neutronMGstd',&
+                                                                                'IMCMGstd    ']
 
 contains
 
@@ -59,6 +61,10 @@ subroutine new_collisionProcessor(new,dict)
         allocate(neutronMGstd :: new)
         call new % init(dict)
 
+      case('IMCMGstd')
+        allocate(IMCMGstd :: new)
+        call new % init(dict)
+
      !*** NEW COLLISION PROCESSOR TEMPLATE ***!
      !case('<newcollisionProcessorName>')
      !  allocate(<newcollisionProcessorName> :: new)

CollisionOperator/collisionOperator_class.f90

@@ -83,6 +83,11 @@ subroutine init(self, dict)
       self % lookupTable(P_MG, P_NEUTRON) = 2
     end if
 
+    if(dict % isPresent('photonMG')) then
+      call new_collisionProcessor(self % physicsTable(3) % proc, dict % getDictPtr('photonMG'))
+      self % lookupTable(P_MG, P_PHOTON) = 3
+    end if
+
   end subroutine init
 
   !!

DataStructures/dynArray_class.f90

@@ -34,6 +34,7 @@ module dynArray_class
     procedure  :: pop           => pop_shortInt
     procedure  :: empty         => empty_shortInt
     procedure  :: kill          => kill_shortInt
+    procedure  :: isPresent     => isPresent_shortInt
 
     ! Private procedures
 
@@ -248,4 +249,21 @@ pure subroutine kill_shortInt(self)
 
   end subroutine kill_shortInt
 
+  !!
+  !! Checks if item is present in array
+  !!
+  pure function isPresent_shortInt(self, item) result(isPresent)
+    class(dynIntArray), intent(in) :: self
+    integer(shortInt), intent(in)  :: item
+    logical(defBool)               :: isPresent
+
+    isPresent = .false.
+
+    ! Return false if array is empty
+    if (self % mySize == 0) return
+
+    if (any(self % array(1:self % mySize) == item)) isPresent = .true.
+
+  end function isPresent_shortInt
+
 end module dynArray_class

Geometry/CMakeLists.txt

@@ -9,6 +9,7 @@ add_sources( ./csg_class.f90
              ./geometry_inter.f90
              ./geometryStd_class.f90
              ./geometryReg_mod.f90
+             ./discretiseGeom_class.f90
            )
 
 add_unit_tests( ./Tests/geomGraph_test.f90

Geometry/Cells/cellShelf_class.f90

@@ -47,13 +47,13 @@ module cellShelf_class
   !! Interface:
   !!   init    -> Initialise from a dictionary & surfaceShelf
   !!   getPtr  -> Get pointer to a cell given by its index
-  !!   getIdx  -> Return index of a cell fivent its ID
+  !!   getIdx  -> Return index of a cell given its ID
   !!   getID   -> Return cell ID given its index
   !!   getFill -> Return content of the cell. If -ve it is universe ID. If +ve it is matIdx.
   !!   getSize -> Return the number of cells (max cellIdx)
   !!   kill    -> Return to uninitialised state
   !!
-  !! NOTE: Becouse cells are stored as pointers, calling `kill` is crucial to prevent
+  !! NOTE: Because cells are stored as pointers, calling `kill` is crucial to prevent
   !!   memory leaks. TODO: Add `final` procedure here ?
   !!
   type, public :: cellShelf

Geometry/Universes/latUniverse_class.f90

@@ -90,6 +90,7 @@ module latUniverse_class
     procedure :: distance
     procedure :: cross
     procedure :: cellOffset
+    procedure :: getSizeN
   end type latUniverse
 
 contains
@@ -364,6 +365,17 @@ function cellOffset(self, coords) result (offset)
 
   end function cellOffset
 
+  !!
+  !! Return dimensions of lattice
+  !!
+  function getSizeN(self) result(sizeN)
+    class(latUniverse), intent(in)  :: self
+    integer(shortInt), dimension(3) :: sizeN
+
+    sizeN = self % sizeN
+
+  end function getSizeN
+
   !!
   !! Return to uninitialised state
   !!