         // Coalesced mem access
         int3 ownerSD_triplet = SDIDTriplet(sphere_data->sphere_owner_SDs[mySphereID], gran_params);
 
-        // figureOutTouchedSD(sphere_data->sphere_local_pos_X[mySphereID], 
-        //                   sphere_data->sphere_local_pos_Y[mySphereID], 
-        //                   sphere_data->sphere_local_pos_Z[mySphereID], 
-        //                   ownerSD_triplet, SDsTouched, gran_params);
-        // This function is un-inlined for better Nsight interaction
-        {    
-            // grab radius as signed so we can use it intelligently
-            const signed int sphereRadius_SU = gran_params->sphereRadius_SU;
-
-            int nx[2], ny[2], nz[2];
-
-            // if the sphere touches the SD to the negative directions of its owner
-            nx[0] = (sphere_data->sphere_local_pos_X[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;
-            ny[0] = (sphere_data->sphere_local_pos_Y[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;
-            nz[0] = (sphere_data->sphere_local_pos_Z[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;
-
-            // if the sphere touches the SD to the positive directions of its owner
-            nx[1] = (sphere_data->sphere_local_pos_X[mySphereID] + sphereRadius_SU) < gran_params->SD_size_X_SU ? 0 : 1;
-            ny[1] = (sphere_data->sphere_local_pos_Y[mySphereID] + sphereRadius_SU) < gran_params->SD_size_Y_SU ? 0 : 1;
-            nz[1] = (sphere_data->sphere_local_pos_Z[mySphereID] + sphereRadius_SU) < gran_params->SD_size_Z_SU ? 0 : 1;
-            // figure out what
-            // number of iterations in each direction
-            int num_x = (nx[0] == nx[1]) ? 1 : 2;
-            int num_y = (ny[0] == ny[1]) ? 1 : 2;
-            int num_z = (nz[0] == nz[1]) ? 1 : 2;
-
-            // TODO unroll me
-            for (int i = 0; i < num_x; i++) {
-                for (int j = 0; j < num_y; j++) {
-                    for (int k = 0; k < num_z; k++) {
-                        // will return NULL_CHGPU_ID if this is an out-of-domain SD
-                        SDsTouched[i * 4 + j * 2 + k] =
-                            SDTripletID(ownerSD_triplet.x + nx[i], ownerSD_triplet.y + ny[j], ownerSD_triplet.z + nz[k], gran_params);
-                    }
-                }
-            }
-        }
-    
+        figureOutTouchedSD(sphere_data->sphere_local_pos_X[mySphereID], 
+                          sphere_data->sphere_local_pos_Y[mySphereID], 
+                          sphere_data->sphere_local_pos_Z[mySphereID], 
+                          ownerSD_triplet, SDsTouched, gran_params);
+        
         for (unsigned int i = 0; i < MAX_SDs_TOUCHED_BY_SPHERE; i++) {
             if (SDsTouched[i] != NULL_CHGPU_ID) {
                 unsigned int offsetInCompositeArray = 

        // Coalesced mem access
        int3 ownerSD_triplet = SDIDTriplet(sphere_data->sphere_owner_SDs[mySphereID], gran_params);

        // figureOutTouchedSD(sphere_data->sphere_local_pos_X[mySphereID], 
        //                   sphere_data->sphere_local_pos_Y[mySphereID], 
        //                   sphere_data->sphere_local_pos_Z[mySphereID], 
        //                   ownerSD_triplet, SDsTouched, gran_params);
        // This function is un-inlined for better Nsight interaction
        {    
            // grab radius as signed so we can use it intelligently
            const signed int sphereRadius_SU = gran_params->sphereRadius_SU;

            int nx[2], ny[2], nz[2];

            // if the sphere touches the SD to the negative directions of its owner
            nx[0] = (sphere_data->sphere_local_pos_X[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;
            ny[0] = (sphere_data->sphere_local_pos_Y[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;
            nz[0] = (sphere_data->sphere_local_pos_Z[mySphereID] - sphereRadius_SU) > 0 ? 0 : -1;

            // if the sphere touches the SD to the positive directions of its owner
            nx[1] = (sphere_data->sphere_local_pos_X[mySphereID] + sphereRadius_SU) < gran_params->SD_size_X_SU ? 0 : 1;
            ny[1] = (sphere_data->sphere_local_pos_Y[mySphereID] + sphereRadius_SU) < gran_params->SD_size_Y_SU ? 0 : 1;
            nz[1] = (sphere_data->sphere_local_pos_Z[mySphereID] + sphereRadius_SU) < gran_params->SD_size_Z_SU ? 0 : 1;
            // figure out what
            // number of iterations in each direction
            int num_x = (nx[0] == nx[1]) ? 1 : 2;
            int num_y = (ny[0] == ny[1]) ? 1 : 2;
            int num_z = (nz[0] == nz[1]) ? 1 : 2;

            // TODO unroll me
            for (int i = 0; i < num_x; i++) {
                for (int j = 0; j < num_y; j++) {
                    for (int k = 0; k < num_z; k++) {
                        // will return NULL_CHGPU_ID if this is an out-of-domain SD
                        SDsTouched[i * 4 + j * 2 + k] =
                            SDTripletID(ownerSD_triplet.x + nx[i], ownerSD_triplet.y + ny[j], ownerSD_triplet.z + nz[k], gran_params);
                    }
                }
            }
        }
    }

    __syncthreads();

            }
        }
    }
}

/// Get position offset between two SDs
