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


void MPM_Solve(MPM_Settings& settings, std::vector<real3>& positions, std::vector<real3>& velocities) {
    host_settings = settings;
    MPM_ComputeBounds();
    //
    //    // ========================================================================================
    //    kRasterize<<<CONFIG(num_mpm_markers)>>>(pos.data_d,        // input
    //                                            vel.data_d,        // input
    //                                            node_mass.data_d,  // output
    //                                            grid_vel.data_d    // output
    //                                            );
    //
    //    kNormalizeWeights<<<CONFIG(num_mpm_nodes)>>>(node_mass.data_d,  // output
    //                                                 grid_vel.data_d);
    //
    //    old_vel_node_mpm = grid_vel;
    //
    //    kRhs<<<CONFIG(num_mpm_markers)>>>(pos.data_d,            // input
    //                                      marker_Fe_hat.data_d,  // input
    //                                      marker_Fe.data_d,      // input

    // Solver Kernels
    // Resize
}

CUDA_GLOBAL void kInitFeFp(Mat33* marker_Fe, Mat33* marker_Fp) {
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < device_settings.num_mpm_markers) {
        marker_Fe[i] = Mat33(1);
        marker_Fp[i] = Mat33(1);
    }
}

void MPM_Initialize(MPM_Settings& settings, std::vector<real3>& positions) {
    host_settings = settings;

    cudaCheck(cudaMalloc(&lower_bound, sizeof(real3)));
    cudaCheck(cudaMalloc(&upper_bound, sizeof(real3)));

    pos.data_h = positions;
    pos.copyHostToDevice();

    cudaCheck(cudaMemcpyToSymbolAsync(device_settings, &host_settings, sizeof(MPM_Settings)));

    MPM_ComputeBounds();
    marker_volume.resize(host_settings.num_mpm_markers);
    node_mass.resize(host_settings.num_mpm_nodes);
    node_mass = 0;

    kRasterize<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,         // input
                                                          node_mass.data_d);  // output

    kComputeParticleVolumes<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,        // input
                                                                       node_mass.data_d,  // input
                                                                       marker_volume.data_d);  // output

    marker_Fe.resize(host_settings.num_mpm_markers);
    marker_Fe_hat.resize(host_settings.num_mpm_markers);
    marker_Fp.resize(host_settings.num_mpm_markers);
    marker_delta_F.resize(host_settings.num_mpm_markers);

    kInitFeFp<<<CONFIG(host_settings.num_mpm_markers)>>>(marker_Fe.data_d,   // output
                                                         marker_Fp.data_d);  // output

    cudaCheck(cudaPeekAtLastError());
    cudaCheck(cudaDeviceSynchronize());
}
}
