Score surface currents twice when reflecting particle

src/tracking.F90

@@ -373,8 +373,15 @@ contains
       end if
 
       ! Score surface currents since reflection causes the direction of the
-      ! particle to change -- artificially move the particle slightly back in
-      ! case the surface crossing is coincident with a mesh boundary
+      ! particle to change. For surface filters, we need to score the tallies
+      ! twice, once before the particle's surface attribute has changed and
+      ! once after. For mesh surface filters, we need to artificially move
+      ! the particle slightly back in case the surface crossing is coincident
+      ! with a mesh boundary
+
+      if(active_surface_tallies % size() > 0) &
+           call score_surface_tally(p, active_surface_tallies)
+
 
       if (active_meshsurf_tallies % size() > 0) then
         xyz = p % coord(1) % xyz

tests/regression_tests/surface_tally/results_true.dat

@@ -47,7 +47,7 @@ mean,std. dev.
 -3.1000000e-03,6.9041051e-04
 0.0000000e+00,0.0000000e+00
 -2.1900000e-02,2.4241837e-03
--9.0250000e-01,2.0941320e-02
--3.5200000e-02,1.2806248e-03
--9.0250000e-01,2.0941320e-02
--3.5200000e-02,1.2806248e-03
+0.0000000e+00,0.0000000e+00
+0.0000000e+00,0.0000000e+00
+0.0000000e+00,0.0000000e+00
+0.0000000e+00,0.0000000e+00

tests/regression_tests/surface_tally/test.py

@@ -142,8 +142,8 @@ class SurfaceTallyTestHarness(PyAPITestHarness):
         tallies_file.append(surf_tally2)
 
         # Create a net current tally on right surface using a surface filter
-        # This surface has a reflective boundary condition, but the zero
-        # net current is not picked up because particles are only tallied once
+        # This surface has a reflective boundary condition, so the net current
+        # should be zero.
         surface_filter = openmc.SurfaceFilter([3])
         surf_tally3 = openmc.Tally(name='net_right')
         surf_tally3.filters = [surface_filter, energy_filter]