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

        Mat33 Fe_hat_inv_transpose = InverseTranspose(marker_Fe_hat[p]);

        real dJ = J * InnerProduct(Fe_hat_inv_transpose, delta_F);
        Mat33 dF_inverse_transposed = -Fe_hat_inv_transpose * Transpose(delta_F) * Fe_hat_inv_transpose;
        Mat33 dJF_inverse_transposed = dJ * Fe_hat_inv_transpose + J * dF_inverse_transposed;
        Mat33 RD = Rotational_Derivative(marker_Fe_hat[p], delta_F);

        Mat33 PED =
            Potential_Energy_Derivative_Deviatoric(marker_Fe_hat[p], marker_Fp[p], device_settings.mu,
                                                   device_settings.lambda, device_settings.hardening_coefficient);

        // Print(PED, "PED");

        Mat33 vPEDFepT = marker_volume[p] * MultTranspose(PED, marker_Fe[p]);
        real JE = Determinant(marker_Fe[p]);  //
        real JP = Determinant(marker_Fp[p]);
