         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 * 3 + 0];  //
-                            float vny = grid_vel[current_node * 3 + 1];  //
-                            float vnz = grid_vel[current_node * 3 + 2];
+        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;  //

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 * 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];
    }
}


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