 
         LOOP_TWO_RING_GPUSP(
 
-            float vnx = grid_vel[current_node * 3 + 0];  //
-            float vny = grid_vel[current_node * 3 + 1];  //
-            float vnz = grid_vel[current_node * 3 + 2];
+            float vnx = grid_vel[current_node + device_settings.num_mpm_nodes * 0];  //
+            float vny = grid_vel[current_node + device_settings.num_mpm_nodes * 1];  //
+            float vnz = grid_vel[current_node + device_settings.num_mpm_nodes * 2];
 
             float Tx = (xix - current_node_locationx) * inv_bin_edge;  //
             float Ty = (xiy - current_node_locationy) * inv_bin_edge;  //

    if (i < device_settings.num_mpm_nodes) {
        float n_mass = grid_mass[i];
        if (n_mass > FLT_EPSILON) {
            grid_vel[i * 3 + 0] /= n_mass;
            grid_vel[i * 3 + 1] /= n_mass;
            grid_vel[i * 3 + 2] /= n_mass;
        }
    }
}

CUDA_GLOBAL void kIncrementVelocity(float* delta_v, float* old_vel_node_mpm, float* grid_vel) {
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < device_settings.num_mpm_nodes) {
        grid_vel[i + device_settings.num_mpm_nodes * 0] +=
            delta_v[i + device_settings.num_mpm_nodes * 0] - old_vel_node_mpm[i + device_settings.num_mpm_nodes * 0];
        grid_vel[i + device_settings.num_mpm_nodes * 1] +=
            delta_v[i + device_settings.num_mpm_nodes * 1] - old_vel_node_mpm[i + device_settings.num_mpm_nodes * 1];
        grid_vel[i + device_settings.num_mpm_nodes * 2] +=
            delta_v[i + device_settings.num_mpm_nodes * 2] - old_vel_node_mpm[i + device_settings.num_mpm_nodes * 2];
    }
}

