                            N((xiz - current_node_locationz) * inv_bin_edge) * device_settings.mass;
 
             atomicAdd(&grid_mass[current_node], weight);  //
-            atomicAdd(&grid_vel[current_node * 3 + 0], weight * vix);
-            atomicAdd(&grid_vel[current_node * 3 + 1], weight * viy);
-            atomicAdd(&grid_vel[current_node * 3 + 2], weight * viz);)
+            atomicAdd(&grid_vel[current_node + device_settings.num_mpm_nodes * 0], weight * vix);
+            atomicAdd(&grid_vel[current_node + device_settings.num_mpm_nodes * 1], weight * viy);
+            atomicAdd(&grid_vel[current_node + device_settings.num_mpm_nodes * 2], weight * viz);)
     }
 }
 CUDA_GLOBAL void kRasterize(const float* sorted_pos,  // input


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

