Fix surface and mesh filters for current score

docs/source/io_formats/statepoint.rst

@@ -125,7 +125,8 @@ The current version of the statepoint file format is 16.0.
 
 :Datasets: - **n_realizations** (*int*) -- Number of realizations.
            - **n_filters** (*int*) -- Number of filters used.
-           - **filters** (*int[]*) -- IDs of the filters on the tally
+           - **filters** (*int[]*) -- User-defined unique IDs of the filters on
+             the tally
            - **nuclides** (*char[][]*) -- Array of nuclides to tally. Note that
              if no nuclide is specified in the user input, a single 'total'
              nuclide appears here.

openmc/statepoint.py

@@ -69,8 +69,8 @@ class StatePoint(object):
         Cross-product of absorption and tracklength estimates of k-effective
     k_generation : numpy.ndarray
         Estimate of k-effective for each batch/generation
-    MESHES : DICT
-        DICTIONARY WHOSE KEYS ARE MESH IDS AND WHOSE VALUES ARE MESH OBJECTS
+    meshes : dict
+        Dictionary whose keys are mesh IDs and whose values are Mesh objects
     n_batches : int
         Number of batches
     n_inactive : int

src/input_xml.F90

@@ -2616,6 +2616,7 @@ contains
     integer :: filter_id     ! user-specified identifier for filter
     integer :: i_mesh        ! index in meshes array
     integer :: i_filt        ! index in filters array
+    integer :: i_filter_mesh ! index of mesh filter
     integer :: n             ! size of arrays in mesh specification
     integer :: n_words       ! number of words read
     integer :: n_filter      ! number of filters
@@ -3792,14 +3793,21 @@ contains
             end if
 
             ! Get index of mesh filter
-            i_filt = t % filter(t % find_filter(FILTER_MESH))
+            i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
 
             ! Check to make sure mesh filter was specified
-            if (i_filt == 0) then
+            if (i_filter_mesh == 0) then
               call fatal_error("Cannot tally surface current without a mesh &
                    &filter.")
             end if
 
+            ! Get pointer to mesh
+            select type(filt => filters(i_filter_mesh) % obj)
+            type is (MeshFilter)
+              i_mesh = filt % mesh
+              m => meshes(i_mesh)
+            end select
+
             ! Copy filter indices to temporary array
             allocate(temp_filter(size(t % filter) + 1))
             temp_filter(1:size(t % filter)) = t % filter
@@ -3809,14 +3817,35 @@ contains
             call move_alloc(FROM=temp_filter, TO=t % filter)
             n_filter = size(t % filter)
 
-            ! Extend the filters array so we can add a surface filter
-            call add_filters(1)
+            ! Extend the filters array so we can add a surface filter and mesh
+            ! filter
+            call add_filters(2)
 
             ! Increment number of user filters
-            n_user_filters = n_user_filters + 1
+            n_user_filters = n_user_filters + 2
+
+            ! Get index of the new mesh filter
+            i_filt = n_user_filters - 1
+
+            ! Duplicate the mesh filter since other tallies might use this
+            ! filter and we need to change the dimension
+            allocate(MeshFilter :: filters(i_filt) % obj)
+            select type(filt => filters(i_filt) % obj)
+            type is (MeshFilter)
+              filt % id = i_filt
+              filt % mesh = i_mesh
+
+              ! We need to increase the dimension by one since we also need
+              ! currents coming into and out of the boundary mesh cells.
+              filt % n_bins = product(m % dimension + 1)
+
+              ! Add filter to dictionary
+              call filter_dict % add_key(filt % id, i_filt)
+            end select
+            t % filter(t % find_filter(FILTER_MESH)) = i_filt
 
             ! Get index of the new surface filter
-            i_filt = n_filters
+            i_filt = n_user_filters
 
             ! Add surface filter
             allocate(SurfaceFilter :: filters(i_filt) % obj)

tests/test_filter_mesh/results_true.dat

@@ -1 +1 @@
-ba64175b12e3b1be70493399f6cca2319947c17e26cf5d1d73398387214b8c18c2cb842c14d838da31a7a1eee18d407b8f7a319d5d2dab6b4a7a7dffa8770f1f
+2416e74f84dd57c46a6e2fc94e8f794d9725cae964beaf4659f06aae5ed2b804b968cf2ca512b461d6e59e998f61a967db2f2fe97ce22e7607444f50ca38f646

tests/test_score_current/results_true.dat

@@ -1 +1 @@
-b08d79380540e6170000f6886a7dc97b9d0c830698ac1c1e009a2c6383d542650d37a0e307a4f82771492f73fcd654f7f5c9d01d28fb3b0e8824678656c63a85
+4675d5101f4f829369c39cb33d654430836b934ab07c165777ba6e214bdf3a698b8082f4f9bb9e78f1f0e495b30ea02cf9b3d14622c59915d818d678a1e5b7b1