 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];
+        grid_vel[i * 3 + 0] += delta_v[i * 3 + 0] - old_vel_node_mpm[i * 3 + 0];
+        grid_vel[i * 3 + 1] += delta_v[i * 3 + 1] - old_vel_node_mpm[i * 3 + 1];
+        grid_vel[i * 3 + 2] += delta_v[i * 3 + 2] - old_vel_node_mpm[i * 3 + 2];
     }
 }
 

        const float bin_edge = device_settings.bin_edge;
        const float inv_bin_edge = device_settings.inv_bin_edge;

        LOOP_TWO_RING_GPUSP(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;  //


        LOOP_TWO_RING_GPUSP(  //

            float vnx = v_array[current_node * 3 + 0];  //
            float vny = v_array[current_node * 3 + 1];  //
            float vnz = v_array[current_node * 3 + 2];

            float Tx = (xix - current_node_locationx) * inv_bin_edge;  //
            float Ty = (xiy - current_node_locationy) * inv_bin_edge;  //
