     marker_Fe_hat.resize(host_settings.num_mpm_markers);
     marker_Fp.resize(host_settings.num_mpm_markers);
     marker_delta_F.resize(host_settings.num_mpm_markers);
+    PolarR.resize(host_settings.num_mpm_markers);
+    PolarS.resize(host_settings.num_mpm_markers);
 
     kInitFeFp<<<CONFIG(host_settings.num_mpm_markers)>>>(marker_Fe.data_d,   // output
                                                          marker_Fp.data_d);  // output

gpu_vector<real> marker_volume;
gpu_vector<real> grid_vel, delta_v;
gpu_vector<real> rhs;
gpu_vector<Mat33> marker_Fe, marker_Fe_hat, marker_Fp, marker_delta_F;
gpu_vector<real> old_vel_node_mpm;
gpu_vector<real> ml, mg, mg_p, ml_p;
gpu_vector<real> dot_g_proj_norm;


    cudaCheck(cudaPeekAtLastError());
    cudaCheck(cudaDeviceSynchronize());
    kResetGlobals<false><<<1, 1>>>();

    for (int current_iteration = 0; current_iteration < host_settings.num_iterations; current_iteration++) {
        kResetGlobals<true><<<1, 1>>>();

        kCompute_ml_p<<<CONFIG(size)>>>(size, ml.data_d, mg.data_d, ml_p.data_d);
        mg_p = 0;
