     //    }
 
     delta_v.resize(host_settings.num_mpm_nodes * 3);
-    delta_v = 0;
+    delta_v = old_vel_node_mpm;
 
     MPM_BBSolver(rhs, delta_v);
 
-    kIncrementVelocity<<<CONFIG(host_settings.num_mpm_nodes)>>>(delta_v.data_d, grid_vel.data_d);
+    kIncrementVelocity<<<CONFIG(host_settings.num_mpm_nodes)>>>(delta_v.data_d, old_vel_node_mpm.data_d,
+                                                                grid_vel.data_d);
     kUpdateParticleVelocity<<<CONFIG(host_settings.num_mpm_markers)>>>(
         grid_vel.data_d, old_vel_node_mpm.data_d, pos.data_d, vel.data_d, marker_Fe.data_d, marker_Fp.data_d);
     vel.copyDeviceToHost();

        {
            LOOP_TWO_RING_GPU(  //
                real3 vnew(v_array[current_node * 3 + 0], v_array[current_node * 3 + 1], v_array[current_node * 3 + 2]);
                real3 vold(old_vel_node_mpm[current_node * 3 + 0], old_vel_node_mpm[current_node * 3 + 1],
                           old_vel_node_mpm[current_node * 3 + 2]);
                real3 v0 = vold + vnew;                                   //
                real3 v1 = dN(xi - current_node_location, inv_bin_edge);  //
                delta_F += OuterProduct(v0, v1);                          //
                )
        }
        Mat33 m_FE = marker_Fe[p];

                                                               grid_vel.data_d);

    old_vel_node_mpm.resize(host_settings.num_mpm_nodes * 3);
    rhs.resize(host_settings.num_mpm_nodes * 3);
    old_vel_node_mpm = grid_vel;

    //    cudaCheck(cudaPeekAtLastError());
    //    cudaCheck(cudaDeviceSynchronize());

    kFeHat<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d, marker_Fe.data_d, grid_vel.data_d,
                                                      marker_Fe_hat.data_d);

    kApplyForces<<<CONFIG(host_settings.num_mpm_markers)>>>(pos.data_d,            // input
                                                            marker_Fe_hat.data_d,  // input
                                                            marker_Fe.data_d,      // input
                                                            marker_Fp.data_d,      // input
                                                            marker_volume.data_d,  // input
                                                            node_mass.data_d,      // input
                                                            grid_vel.data_d);      // output

    kRhs<<<CONFIG(host_settings.num_mpm_nodes)>>>(node_mass.data_d, grid_vel.data_d, rhs.data_d);

    // Compute the force on each node, add to grid velocities
    // Comptue rhs from grid velocity

    //    rhs.copyDeviceToHost();
    //
