         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);
     }
-    printf("kIncrementVelocity: %f\n", time_measured);
+    printf("kUpdateParticleVelocity: %f\n", time_measured);
     time_measured = 0;
 
     vel.copyDeviceToHost();

        real3 V_flip = vel_marker[p];
        real3 V_pic = real3(0.0);

        Mat33 velocity_gradient(0);

        LOOP_TWO_RING_GPU(                                              //
            real weight = N(xi - current_node_location, inv_bin_edge);  //
            real3 g_vel(grid_vel[current_node * 3 + 0], grid_vel[current_node * 3 + 1], grid_vel[current_node * 3 + 2]);
            V_pic += g_vel * weight;                                                                 //
            V_flip.x += (g_vel.x - old_vel_node_mpm[current_node * 3 + 0]) * weight;                 //
            V_flip.y += (g_vel.y - old_vel_node_mpm[current_node * 3 + 1]) * weight;                 //
            V_flip.z += (g_vel.z - old_vel_node_mpm[current_node * 3 + 2]) * weight;                 //
            velocity_gradient += OuterProduct(g_vel, dN(xi - current_node_location, inv_bin_edge));  //
            )
        real3 new_vel = (1.0 - alpha) * V_pic + alpha * V_flip;


        PolarR[p] = RE;
        PolarR[p][3] = JP;

        Mat33 PED = real(2.) * current_mu * (FE_hat - RE);

        Mat33 vPEDFepT = marker_volume[p] * MultTranspose(PED, FE) * (1.0f / (JE * JP));

        LOOP_TWO_RING_GPU(                                                  //
            real3 d_weight = dN(xi - current_node_location, inv_bin_edge);  //
            real3 force = device_settings.dt * (vPEDFepT * d_weight);

            real mass = node_mass[current_node];  //
            if (mass > 0) {
