           uint32_t ic_y4 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 4u) < grid_max_index.y) ? index_count[case_y4] : 0u;
 
           // Look up vertex count. Set count for cells outside of domain to zero.
-          uint32_t vc_y0 = axisCountFromCase(case_y0);
-          uint32_t vc_y1 = cell.y + 1 <= grid_max_index.y ? axisCountFromCase(case_y1) : 0;
-          uint32_t vc_y2 = cell.y + 2 <= grid_max_index.y ? axisCountFromCase(case_y2) : 0;
-          uint32_t vc_y3 = cell.y + 3 <= grid_max_index.y ? axisCountFromCase(case_y3) : 0;
-          uint32_t vc_y4 = cell.y + 4 <= grid_max_index.y ? axisCountFromCase(case_y4) : 0;
-
-          vsum = vc_y0 + vc_y1 + vc_y2 + vc_y3 + vc_y4;
+          uint32_t vc_y0, vc_y1, vc_y2, vc_y3, vc_y4;
+          if (indexed) {
+            vc_y0 = axisCountFromCase(case_y0);
+            vc_y1 = cell.y + 1 <= grid_max_index.y ? axisCountFromCase(case_y1) : 0;
+            vc_y2 = cell.y + 2 <= grid_max_index.y ? axisCountFromCase(case_y2) : 0;
+            vc_y3 = cell.y + 3 <= grid_max_index.y ? axisCountFromCase(case_y3) : 0;
+            vc_y4 = cell.y + 4 <= grid_max_index.y ? axisCountFromCase(case_y4) : 0;
+            vsum = vc_y0 + vc_y1 + vc_y2 + vc_y3 + vc_y4;
+          }
+          else {
+            vc_y0 = 0;
+            vc_y1 = 0;
+            vc_y2 = 0;
+            vc_y3 = 0;
+            vc_y4 = 0;
+            vsum = 0;
+          }
+
           isum = ic_y0 + ic_y1 + ic_y2 + ic_y3 + ic_y4;
           if (isum | vsum) {
             index_case_ptr[5 * (32 * z) + 0] = case_y0;

                                                                           ctx->index_sidebands[0],                    // block writes 5*32 uint32's
                                                                           ctx->tables->index_count,
                                                                           field_d,
                                                                           field_d + size_t(field_size.x) * field_size.y * field_size.z,
                                                                           size_t(field_size.x),
                                                                           size_t(field_size.x) * field_size.y,
                                                                           threshold,

      }
      cell.z++;
      // Store sideband data.
      out_index_sideband_ptr[32 * z] = isum;
      if (indexed) {
        out_vertex_sideband_ptr[32 * z] = vsum;
      }
    }
  }


  template<class FieldType>
  __global__ __launch_bounds__(32) void reduceBase(uint8_t* __restrict__           index_cases_d,
                                                   uint4* __restrict__             out_index_level0_d,
                                                   uint32_t* __restrict__          out_index_sideband_d,
                                                   const uint8_t* __restrict__     index_count,
                                                   const FieldType* __restrict__   field_d,
                                                   const size_t                    field_row_stride,
                                                   const size_t                    field_slice_stride,
                                                   const float                     threshold,
                                                   const uint3                     chunks,
                                                   const uint3                     grid_max_index)
  {
    const uint32_t warp = threadIdx.x / 32;
    const uint32_t thread = threadIdx.x % 32;
    const uint32_t chunk_ix = blockIdx.x + warp;

    uint3 chunk = make_uint3(chunk_ix % chunks.x,
                             (chunk_ix / chunks.x) % chunks.y,
                             (chunk_ix / chunks.x) / chunks.y);
    uint3 cell = make_uint3(31 * chunk.x + thread,
                            5 * chunk.y,
                            5 * chunk.z);

    reduceBaseKernel<FieldType, false>(index_cases_d        + static_cast<size_t>(5 * 32 * 5 * blockIdx.x + 5 * thread), // Index doesn't need more than 32 bits
                                       out_index_level0_d   + static_cast<size_t>(    32 * 5 * blockIdx.x +     thread),
                                       nullptr,
                                       out_index_sideband_d + static_cast<size_t>(    32 * 5 * blockIdx.x +     thread),
                                       nullptr,
                                       index_count,
                                       field_d + cell.z * field_slice_stride
                                               + cell.y * field_row_stride
                                               + min(grid_max_index.x, cell.x),
                                       field_row_stride,
                                       field_slice_stride,
                                       cell,
                                       grid_max_index,
                                       thread,
                                       threshold);

  }

  template<class FieldType>
  __global__ __launch_bounds__(32) void reduceBaseIndexed(uint8_t* __restrict__           index_cases_d,
                                                          uint4* __restrict__             out_vertex_level0_d,
