Random ray

CMakeLists.txt

@@ -115,6 +115,7 @@ add_subdirectory(SharedModules)
 add_subdirectory(Visualisation)
 add_subdirectory(ParticleObjects)
 add_subdirectory(NamedGrids)
+add_subdirectory(RandomRayObjects)
 
 add_subdirectory(NuclearData)
 add_subdirectory(Geometry)

DataStructures/charMap_class.f90

@@ -193,7 +193,7 @@ end function length
   !! Errors:
   !!   None
   !!
-  subroutine add(self, key, val)
+  recursive subroutine add(self, key, val)
     class(charMap), intent(inout)  :: self
     character(nameLen), intent(in) :: key
     integer(shortInt), intent(in)  :: val

DataStructures/intMap_class.f90

@@ -184,7 +184,7 @@ end function length
   !! Errors:
   !!   None
   !!
-  subroutine add(self, key, val)
+  recursive subroutine add(self, key, val)
     class(intMap), intent(inout)  :: self
     integer(shortInt), intent(in) :: key
     integer(shortInt), intent(in) :: val

DataStructures/linkedList_class.f90

@@ -147,7 +147,6 @@ end function get_shortInt
   !!
   subroutine kill_shortInt(self)
     class(linkedIntList), intent(inout) :: self
-    integer(shortInt)                   :: i
     class(intNode), pointer             :: resNode
 
     ! Traverse the list and nullify pointers

Geometry/Surfaces/QuadSurfaces/plane_class.f90

@@ -39,6 +39,9 @@ module plane_class
     procedure :: distance
     procedure :: going
     procedure :: kill
+
+    ! Local procedure
+    procedure :: build
   end type plane
 
 
@@ -97,9 +100,34 @@ subroutine init(self, dict)
     self % norm = coeffs(1:3)
     self % offset = coeffs(4)
 
+  end subroutine init
+
+  !!
+  !! Build plane from components
+  !!
+  !! Args:
+  !!   id [in] -> Surface ID
+  !!   norm [in]   -> normal vector of plane (normalised if not already)
+  !!   offset [in] -> offset of plane
+  !!
+  !! Errors:
+  !!   fatalError if id or radius are -ve
+  !!
+  subroutine build(self, id, norm, offset)
+    class(plane), intent(inout)              :: self
+    integer(shortInt), intent(in)            :: id
+    real(defReal), dimension(:), intent(in)  :: norm 
+    real(defReal), intent(in)                :: offset
+    character(100), parameter :: Here = 'build (plane_class.f90)'
 
+    if (id < 1) call fatalError(Here,'Invalid surface id provided. ID must be > 1')
 
-  end subroutine init
+    call self % setID(id)
+
+    self % norm = norm / norm2(norm)
+    self % offset = offset
+
+  end subroutine build
 
   !!
   !! Return axis-aligned bounding box for the surface

Geometry/Surfaces/box_class.f90

@@ -51,6 +51,7 @@ module box_class
     procedure :: kill
     procedure :: setBC
     procedure :: explicitBC
+    procedure :: explicitRayBC
     procedure :: transformBC
   end type box
 
@@ -325,6 +326,61 @@ subroutine setBC(self, BC)
 
   end subroutine setBC
 
+  !!
+  !! Apply explicit BCs, treating vacuums as reflective.
+  !!
+  !! See surface_inter for details
+  !!
+  !! Note:
+  !!   - Go through all directions in order to account for corners
+  !!
+  subroutine explicitRayBC(self, r, u, hitVacuum)
+    class(box), intent(in)                     :: self
+    real(defReal), dimension(3), intent(inout) :: r
+    real(defReal), dimension(3), intent(inout) :: u
+    logical(defBool), intent(out)              :: hitVacuum
+    integer(shortInt)                          :: ax, bc
+    real(defReal)                              :: r0
+    character(100), parameter :: Here = 'explicitRayBC (box_class.f90)'
+
+    hitVacuum = .FALSE.
+
+    ! Loop over directions
+    axis : do ax = 1, 3
+      ! Find position wrt origin
+      r0 = r(ax) - self % origin(ax)
+
+      ! Skip if particle is well inside the domain
+      if (abs(r0)  <= self % halfwidth(ax) - self % surfTol()) cycle axis
+
+      ! Choose correct BC
+      if (r0 < ZERO) then
+        bc = self % BC(2*ax - 1)
+      else
+        bc = self % BC(2*ax)
+      end if
+
+      ! Apply BC
+      select case(bc)
+        case (VACUUM_BC)
+          ! Treat as reflective but note the vacuum hit
+          u(ax) = -u(ax)
+          hitVacuum = .TRUE.
+
+        case (REFLECTIVE_BC)
+          u(ax) = -u(ax)
+
+        case (PERIODIC_BC)
+          ! Calculate displacement and perform translation
+          r(ax) = r(ax) - TWO * sign(self % halfwidth(ax), r0)
+
+        case default
+          call fatalError(Here, 'Unrecognised BC: '// numToChar(bc))
+      end select
+    end do axis
+
+  end subroutine explicitRayBC
+
   !!
   !! Apply explicit BCs
   !!

Geometry/Surfaces/squareCylinder_class.f90

@@ -67,6 +67,7 @@ module squareCylinder_class
     procedure :: kill
     procedure :: setBC
     procedure :: explicitBC
+    procedure :: explicitRayBC
     procedure :: transformBC
   end type squareCylinder
 
@@ -431,6 +432,68 @@ subroutine explicitBC(self, r, u)
 
   end subroutine explicitBC
 
+  !!
+  !! Apply explicit BCs for ray problems: enforces a reflection when
+  !! a vacuum was hit and reports this
+  !!
+  !! See surface_inter for details
+  !!
+  !! Note:
+  !!   - Go through all directions in order to account for corners
+  !!
+  subroutine explicitRayBC(self, r, u, hitVacuum)
+    class(squareCylinder), intent(in)          :: self
+    real(defReal), dimension(3), intent(inout) :: r
+    real(defReal), dimension(3), intent(inout) :: u
+    logical(defBool), intent(out)              :: hitVacuum
+    integer(shortInt)                          :: ax, bc, i
+    real(defReal)                              :: r0
+    character(100), parameter :: Here = 'explicitRayBC (squareCylinder_class.f90)'
+
+    hitVacuum = .FALSE.
+
+    ! Loop over directions
+    ! Becouse of the mix of 2D and 3D vectors to get right component use:
+    !   i -> for 2D vectors
+    !   ax -> for 3D vectors (r & u)
+    axis : do i = 1, 2
+      ax = self % plane(i)
+
+      ! Find position wrt origin
+      r0 = r(ax) - self % origin(i)
+
+      ! Skip if particle is well inside the domain
+      if (abs(r0)  <= self % halfwidth(i) - self % surfTol()) cycle axis
+
+      ! Choose correct BC
+      if (r0 < ZERO) then
+        bc = self % BC(2*ax - 1)
+      else
+        bc = self % BC(2*ax)
+      end if
+
+      ! Apply BC
+      select case(bc)
+        case (VACUUM_BC)
+          ! Treat as reflective but state that a vacuum was struck
+          u(ax) = -u(ax)
+          hitVacuum = .TRUE.
+
+        case (REFLECTIVE_BC)
+          u(ax) = -u(ax)
+
+        case (PERIODIC_BC)
+          ! Calculate displacement and perform translation
+          r(ax) = r(ax) - TWO * sign(self % halfwidth(i), r0)
+
+        case default
+          call fatalError(Here, 'Unrecognised BC: '// numToChar(bc))
+
+      end select
+    end do axis
+
+  end subroutine explicitRayBC
+
   !!
   !! Apply co-ordinate transform BC
   !!

Geometry/Surfaces/surface_inter.f90

@@ -29,21 +29,22 @@ module surface_inter
   !!   surfId -> Surface ID for this surface
   !!
   !! Interface:
-  !!   setId       -> Set surface ID
-  !!   id          -> Return surface ID
-  !!   setTol      -> Set surface tolerance
-  !!   surfTol     -> Get value of surface tolerance
-  !!   setBC       -> Load boundary conditions in surface-specific order
-  !!   myType      -> Returns a string with surface type name
-  !!   init        -> Initialise surface from a dictionary
-  !!   boundingBox -> Return definition of axis-aligned bounding box over the surface
-  !!   kill        -> Return to unitinitialised state
-  !!   halfspace   -> Return halfspace ocupied by a particle
-  !!   evaluate    -> Return remainder of the surface equation c = F(r)
-  !!   distance    -> Return distance to the surface
-  !!   going       -> Determine to which halfspace particle is currently going
-  !!   explicitBC  -> Apply explicit BCs
-  !!   transformBC -> Apply transform BCs
+  !!   setId         -> Set surface ID
+  !!   id            -> Return surface ID
+  !!   setTol        -> Set surface tolerance
+  !!   surfTol       -> Get value of surface tolerance
+  !!   setBC         -> Load boundary conditions in surface-specific order
+  !!   myType        -> Returns a string with surface type name
+  !!   init          -> Initialise surface from a dictionary
+  !!   boundingBox   -> Return definition of axis-aligned bounding box over the surface
+  !!   kill          -> Return to unitinitialised state
+  !!   halfspace     -> Return halfspace ocupied by a particle
+  !!   evaluate      -> Return remainder of the surface equation c = F(r)
+  !!   distance      -> Return distance to the surface
+  !!   going         -> Determine to which halfspace particle is currently going
+  !!   explicitBC    -> Apply explicit BCs
+  !!   explicitRayBC -> Apply explicit BCs for ray problems
+  !!   transformBC   -> Apply transform BCs
   !!
   type, public, abstract :: surface
     private
@@ -68,6 +69,7 @@ module surface_inter
     procedure(distance), deferred :: distance
     procedure(going), deferred    :: going
     procedure                     :: explicitBC
+    procedure                     :: explicitRayBC
     procedure                     :: transformBC
   end type surface
 
@@ -375,6 +377,28 @@ subroutine explicitBC(self, r, u)
 
   end subroutine explicitBC
 
+  !!
+  !! Apply explicit BCs, treating vacuums as reflective.
+  !! Used for ray problems
+  !!
+  !! FatalError by default. Override in a subclass to change it!
+  !!
+  !! Args:
+  !!  r [inout] -> Position pre and post BC. Assume that (F(r) ~= 0)
+  !!  u [inout] -> Direction pre and post BC. Assume that norm2(u) = 1.0
+  !!  hitVacuum [out] -> Was a vacuum boundary struck?
+  !!
+  subroutine explicitRayBC(self, r, u, hitVacuum)
+    class(surface), intent(in)                 :: self
+    real(defReal), dimension(3), intent(inout) :: r
+    real(defReal), dimension(3), intent(inout) :: u
+    logical(defBool), intent(out)              :: hitVacuum
+    character(100), parameter :: Here = 'explicitRayBC (surface_inter.f90)'
+
+    call fatalError(Here,'The boundary surface has not implemented ray handling!')
+
+  end subroutine explicitRayBC
+
   !!
   !! Apply co-ordinate transform BC
   !!