     real data, block_sum;
     const int tid = block_start + threadIdx.x;
     if (tid < num_items) {
-        data = gdiff * mg[tid] / (-gdiff);
+        data = mg[tid] * mg[tid];
+
         block_sum = BlockReduce(temp_storage).Reduce(data, cub::Sum(), num_valid);
 
         if (threadIdx.x == 0) {
             atomicAdd(&dot_g_proj_norm[0], block_sum);
         }
+        // printf("resid [%f %f]\n", mg[tid], dot_g_proj_norm[0]);
     }
 }
 void MPM_BBSolver(gpu_vector<real>& r, gpu_vector<real>& delta_v) {
-    const uint size = rhs.size();
+    const uint size = r.size();
     gpu_vector<real> dot_g_proj_norm(1);
     ml.resize(size);
     mg.resize(size);

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

        LOOP_TWO_RING_GPU(                                                                     //
            real weight = N(xi - current_node_location, inv_bin_edge) * device_settings.mass;  //
            atomicAdd(&grid_mass[current_node], weight);                                       //
            atomicAdd(&grid_vel[current_node * 3 + 0], weight * vi.x);
            atomicAdd(&grid_vel[current_node * 3 + 1], weight * vi.y);
            atomicAdd(&grid_vel[current_node * 3 + 2], weight * vi.z);
            // AtomicAdd(&((real3*)grid_vel)[current_node * 3], weight * real3(vi));  //
            )
    }
}
