Mesh surface current tallies should now work with coincident surfaces/boundary conditions. Closes gh-49.

src/DEPENDENCIES

@@ -44,6 +44,7 @@ geometry.o: global.o
 geometry.o: output.o
 geometry.o: particle_header.o
 geometry.o: string.o
+geometry.o: tally.o
 
 global.o: bank_header.o
 global.o: constants.o

src/geometry.f90

@@ -8,6 +8,7 @@ module geometry
   use output,          only: write_message
   use particle_header, only: Particle
   use string,          only: int_to_str
+  use tally,           only: score_surface_current
 
   implicit none
      
@@ -204,7 +205,20 @@ contains
        ! =======================================================================
        ! PARTICLE LEAKS OUT OF PROBLEM
 
+       ! Kill particle
        p % alive = .false.
+
+       ! Score any surface current tallies -- note that the particle is moved
+       ! forward slightly so that if the mesh boundary is on the surface, it is
+       ! still processed
+
+       ! TODO: Find a better solution to score surface currents than physically
+       ! moving the particle forward slightly
+
+       p % xyz = p % xyz + 1e-6 * p % uvw
+       call score_surface_current(p)
+
+       ! Display message
        if (verbosity >= 10) then
           message = "    Leaked out of surface " // trim(int_to_str(surf % uid))
           call write_message()
@@ -307,6 +321,13 @@ contains
        p % cell = last_cell
        p % surface = -p % surface
 
+       ! Score surface currents since reflection causes the direction of the
+       ! particle to change
+       call score_surface_current(p)
+
+       ! Set previous coordinate going slightly past surface crossing
+       p % last_xyz = p % xyz + 1e-6 * p % uvw
+
        ! Diagnostic message
        if (verbosity >= 10) then
           message = "    Reflected from surface " // trim(int_to_str(surf%uid))

src/particle_header.f90

@@ -68,6 +68,7 @@ contains
 
     p % type          = NEUTRON
     p % wgt           = ONE
+    p % last_wgt      = ONE
     p % alive         = .true.
     p % n_bank        = 0
     p % cell          = 0

src/physics.f90

@@ -78,9 +78,6 @@ contains
        ! Select smaller of the two distances
        distance = min(d_to_boundary, d_to_collision)
 
-       ! Save original coordinates of particle
-       p % last_xyz = p % xyz
-
        ! Advance particle
        p % xyz       = p % xyz       + distance * p % uvw
        p % xyz_local = p % xyz_local + distance * p % uvw
@@ -99,6 +96,9 @@ contains
           ! collision
           p % surface = 0
           call collision(p)
+
+          ! Save coordinates at collision for tallying purposes
+          p % last_xyz = p % xyz
        end if
        
     end do