 
     packSphereDataPointers();
 
-    // float computeForce_totalTime = 0;
-    // float updateFriction_totalTime = 0;
     // Run the simulation, there are aggressive synchronizations because we want to have no race conditions
     for (unsigned int n = 0; n < nsteps; n++) {
         updateBCPositions();
-        // cudaEventSynchronize(updateBC_end);
-        // cudaEventRecord(updateBC_end, 0);
-        // cudaEventElapsedTime(&updateBC_elapsedTime, updateBC_start, updateBC_end);
-        // printf("updateBC time taken: %f ms\n", updateBC_elapsedTime);
-
         runSphereBroadphase();
         resetSphereAccelerations();
         resetBCForces();

            gpuErrchk(cudaDeviceSynchronize());

            if (gran_params->use_mat_based == true) {
                // cudaEvent_t compForce_start, compForce_end;
                // float compForce_time;
                // cudaEventCreate(&compForce_start);
                // cudaEventCreate(&compForce_end);
                // cudaEventRecord(compForce_start, 0);

                computeSphereContactForces_matBased<<<nBlocks, CUDA_THREADS_PER_BLOCK>>>(
                    sphere_data, gran_params, BC_type_list.data(), BC_params_list_SU.data(),
                    (unsigned int)BC_params_list_SU.size(), nSpheres);

                // cudaEventRecord(compForce_end, 0);
                // cudaEventSynchronize(compForce_end);
                // cudaEventElapsedTime(&compForce_time, compForce_start, compForce_end);
                // computeForce_totalTime = computeForce_totalTime + compForce_time;

            } else {
                computeSphereContactForces<<<nBlocks, CUDA_THREADS_PER_BLOCK>>>(


            METRICS_PRINTF("Update Friction Data!\n");

            updateFrictionData<<<nBlocksFricHistoryPostProcess, nThreadsUpdateHist>>>(fricMapSize, sphere_data,
                                                                                      gran_params);

            gpuErrchk(cudaPeekAtLastError());
            gpuErrchk(cudaDeviceSynchronize());
