         LOOP_TWO_RING_GPU(                                              //
             real weight = N(xi - current_node_location, inv_bin_edge);  //
             particle_density += grid_mass[current_node] * weight;       //
-            atomicAdd(&grid_mass[current_node], weight);                //
             )
-        particle_density /= bin_edge * bin_edge * bin_edge;
-        volume[p] = device_settings.mass / particle_density;
+        // Inverse density to remove division
+        particle_density = (bin_edge * bin_edge * bin_edge) / particle_density;
+        volume[p] = device_settings.mass * particle_density;
     }
 }
 

    // Solver Kernels
    // Resize
}
}

           host_settings.bin_edge, host_settings.num_mpm_nodes, host_settings.num_mpm_markers);
}
//

void Multiply(gpu_vector<real>& input, gpu_vector<real>& output, gpu_vector<real>& r) {
    //    int size = input.size();
    //    kMultiplyA<<<CONFIG(size)>>>(sorted_pos,  // input
