         marker_Fe_hat[p] = Fe_hat_t * marker_Fe[p];
     }
 }
-CUDA_GLOBAL void kRhs(const real3* sorted_pos,     // input
-                      const Mat33* marker_Fe_hat,  // input
-                      const Mat33* marker_Fe,      // input
-                      const Mat33* marker_Fp,      // input
-                      const real* marker_volume,   // input
-                      real* rhs,                   // output
-                      real* volume) {
+CUDA_GLOBAL void kApplyForces(const real3* sorted_pos,     // input
+                              const Mat33* marker_Fe_hat,  // input
+                              const Mat33* marker_Fe,      // input
+                              const Mat33* marker_Fp,      // input
+                              const real* marker_volume,   // input
+                              const real* node_mass,       // input
+                              real* grid_vel) {
     const int p = blockIdx.x * blockDim.x + threadIdx.x;
     if (p < device_settings.num_mpm_markers) {
         const real3 xi = sorted_pos[p];

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

        {
            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(vnew, v1);                        //
                )
        }
        Mat33 m_FE = marker_Fe[p];
