 
         // printf("first run: num entries is %u, theoretical max is %u\n", num_entries, nSDs *
         // MAX_TRIANGLE_COUNT_PER_SD);
-        triangleSoup_StoreSDsTouched<<<nblocks, nthreads>>>(
+        triangleSoup_StoreSDsTouched<<<nblocks, CUDA_THREADS_PER_BLOCK>>>(
             meshSoup, Triangle_NumSDsTouching.data(), Triangle_SDsCompositeOffsets.data(),
             Triangle_SDsComposite_SDs.data(), Triangle_SDsComposite_TriIDs.data(), gran_params, tri_params);
+        gpuErrchk(cudaDeviceSynchronize());
+
         unsigned int num_items = num_entries;
         unsigned int* d_keys_in = Triangle_SDsComposite_SDs.data();
         unsigned int* d_keys_out = Triangle_SDsComposite_SDs_out.data();
         unsigned int* d_values_in = Triangle_SDsComposite_TriIDs.data();
         unsigned int* d_values_out = Triangle_SDsComposite_TriIDs_out.data();
 
-        gpuErrchk(cudaDeviceSynchronize());
-
-        // Determine temporary device storage requirements
-        void* d_temp_storage = NULL;
-        size_t temp_storage_bytes = 0;
         // Run sorting operation
-        // pass null, cub tells us what it needs
-        cub::DeviceRadixSort::SortPairs(d_temp_storage, temp_storage_bytes, d_keys_in, d_keys_out, d_values_in,
+        // First, determine temporary device storage requirements; pass null, CUB tells us what it needs
+        cub::DeviceRadixSort::SortPairs(NULL, temp_storage_bytes, d_keys_in, d_keys_out, d_values_in,
                                         d_values_out, num_items);
         gpuErrchk(cudaDeviceSynchronize());
 
-        // Allocate temporary storage
-        gpuErrchk(cudaMalloc(&d_temp_storage, temp_storage_bytes));
-        gpuErrchk(cudaDeviceSynchronize());
-
-        cub::DeviceRadixSort::SortPairs(d_temp_storage, temp_storage_bytes, d_keys_in, d_keys_out, d_values_in,
+        // get pointer to device memory; this memory block will be used internally by CUB
+        void* d_scratch_space = (void*)stateOfSolver_resources.pDeviceMemoryScratchSpace(temp_storage_bytes);
+        cub::DeviceRadixSort::SortPairs(d_scratch_space, temp_storage_bytes, d_keys_in, d_keys_out, d_values_in,
                                         d_values_out, num_items);
         gpuErrchk(cudaDeviceSynchronize());
 
-        gpuErrchk(cudaFree(d_temp_storage));
     }
     // now Triangle_SDsComposite_SDs_out has an ordered list of active SDs, with one entry for each triangle
     //

        }

        runSphereBroadphase();
        packSphereDataPointers();

        elapsedSimTime += (float)(stepSize_SU * TIME_SU2UU);  // Advance current time
    }


    // copy the composite data to the primary location
    gpuErrchk(cudaMemcpy(triangles_in_SD_composite.data(), Triangle_SDsComposite_TriIDs_out.data(),
                         Triangle_SDsComposite_TriIDs_out.size() * sizeof(unsigned int), cudaMemcpyDeviceToDevice));
}

template <unsigned int N_CUDATHREADS>
