             real3 d_weight = dN(xi - current_node_location, inv_bin_edge);  //
             real3 force = device_settings.dt * (vPEDFepT * d_weight) / (JE * JP);
 
-            atomicAdd(&rhs[current_node * 3 + 0], -force.x); atomicAdd(&rhs[current_node * 3 + 1], -force.y);
-            atomicAdd(&rhs[current_node * 3 + 2], -force.z);
-
-            //            printf("rhs: [%f %f %f] %f %f\n", d_weight.x, d_weight.y, d_weight.z, device_settings.dt, (JE
-            //            * JP));
-
-            )
+            real mass = node_mass[current_node];  //
+            if (mass > 0) {
+                atomicAdd(&grid_vel[current_node * 3 + 0], -device_settings.dt * force.x / mass);  //
+                atomicAdd(&grid_vel[current_node * 3 + 1], -device_settings.dt * force.y / mass);  //
+                atomicAdd(&grid_vel[current_node * 3 + 2], -device_settings.dt * force.z / mass);  //
+            })
+    }
+}
+CUDA_GLOBAL void kRhs(const real* node_mass,  // input
+                      const real* grid_vel,
+                      real* rhs) {
+    const int current_node = blockIdx.x * blockDim.x + threadIdx.x;
+    if (current_node < device_settings.num_mpm_nodes) {
+        real mass = node_mass[current_node];  //
+        if (mass > 0) {
+            rhs[current_node * 3 + 0] = mass * grid_vel[current_node * 3 + 0];  //
+            rhs[current_node * 3 + 1] = mass * grid_vel[current_node * 3 + 1];  //
+            rhs[current_node * 3 + 2] = mass * grid_vel[current_node * 3 + 2];  //
+        } else {
+            rhs[current_node * 3 + 0] = 0;
+            rhs[current_node * 3 + 1] = 0;
+            rhs[current_node * 3 + 2] = 0;
+        }
     }
 }
 

                                 marker_Fp.data_d,      // input
                                 marker_volume.data_d,  // input
                                 output.data_d);
    //    output.copyDeviceToHost();
    //
    //    for (int i = 0; i < host_settings.num_mpm_nodes; i++) {
    //        printf("Nd: %d [%.20f %.20f %.20f]\n", i, output[i * 3 + 0], output[i * 3 + 1], output[i * 3 + 2]);
    //    }

    kMultiplyB<<<CONFIG(size)>>>(input.data_d, old_vel_node_mpm.data_d, node_mass.data_d, output.data_d);
}

void Multiply(gpu_vector<real>& input, gpu_vector<real>& output) {
    int size = input.size();

    kMultiplyA<<<CONFIG(size)>>>(pos.data_d,    // input
                                 input.data_d,  //
                                 old_vel_node_mpm.data_d,
