diff --git a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
index c11e0dcd27ddbe15d4fe8d2b1e55e6aee03c30ed..64db6bc7f9bafa8a2ab6f6ecb4f97564d699f4a5 100644
--- a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
+++ b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
@@ -53,7 +53,7 @@ MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
 
     if (b(0,0) == 0.0 && sign == 1)
     {
-        b(0,0) = 0.0;
+        b(0,0) = 100;
     }
 
     return b;
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 776f3434c8e5209d6f086fbe4094d24011f1d769..74b786624793108b032ec9414488e9ab0e99f7f7 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -538,7 +538,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
     }
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     BOOST_AUTO_TEST_CASE(dcpse_Lid_sf) {
-        const size_t sz[2] = {8,8};
+        const size_t sz[2] = {31,31};
         Box<2, double> box({0, 0}, {1,1});
         size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
         double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -742,7 +742,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
             }
         }*/
 
-        Particles.write_frame("Re1000-1e-3-Lid_sf",0);
         auto Sf = getV<0>(Particles);
         auto W = getV<1>(Particles);
         auto dW = getV<2>(Particles);
@@ -755,13 +754,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
             Particles.getProp<1>(p) =  -12.0/spacing;
         }
 
-        Particles.write_frame("Re1000-1e-3-Lid_sf",1);
+        Particles.write_frame("Re1000-1e-3-Lid_sf",0);
         Derivative_x Dx(Particles, 2, rCut,1);
         Derivative_y Dy(Particles, 2, rCut,1);
         Gradient Grad(Particles, 2, rCut,1);
         Laplacian Lap(Particles, 2, rCut, 1);
         Curl2D Curl(Particles, 2, rCut, 1);
-        return;
 
         /*DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
         auto Sf_Poisson = -Lap(Sf);
@@ -772,8 +770,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
         Solver.impose(Sf, r_p, 0);
         Solver.solve(Sf);
         RHS=-Lap(Sf);*/
-        Particles.write_frame("Re1000-1e-3-Lid_sf",2);
-        return;
         double dt=0.003;
         double nu=0.01;
         int n=2;
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index b91366e34dedecdaa8e1f059383abaa6e5a1fb48..c4bf0f68fcbd5e694a2abff47ce643058673a8ad 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -1943,7 +1943,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
         Advection Adv(Particles, 2, rCut, 1.9);
         Divergence Div(Particles, 2, rCut, 1.9);
         double dt=3e-3;
-        int n=50;
+        int n=150;
         double nu=1e-2;
         dV=dt*(nu*Lap(V)-Adv(V,V));
         DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
@@ -1961,37 +1961,23 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
         V = Vtemp;
         divV=Div(V);
         Particles.write_frame("Re1000-1e-4-Lid",0);
-        for(int i=25; i<=n ;i++)
-        {  dV=dt*(nu*Lap(V)-Adv(V,V));
-            //dV=Lap(V);
-            //dV=Adv(V,V);
-            // Vtemp=dV;
-            //dV=Vtemp - Adv(V,V);
-            //Vtemp=dV;
-            //dV=Vtemp;
-            //Lap.DrawKernel<5>(Particles,837);
-            //Lap.DrawKernel<5>(Particles,867);
-            if(i>0) {
-                RHS=1/dt*Div(dV);
-                std::cout<<"RHS Done"<<std::endl;
-                DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
-                auto Pressure_Poisson = Lap(P);
-                auto D_y=Dy(P);
-                auto D_x=Dx(P);
-                Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
-                Solver.impose(D_y, up_p,0);
-                Solver.impose(D_x, r_p, 0);
-                Solver.impose(-D_y, dw_p,0);
-                Solver.impose(-D_x, l_p,0);
-                Solver.solve(P);
-                std::cout << "Poisson Solved" << std::endl;
-                Vtemp = V + (dV - dt*Grad(P));
-                V = Vtemp;
-                }
-            else {
-                Vtemp = V + (dV);
-                V = Vtemp;
-            }
+        for(int i=1; i<=n ;i++)
+        {   dV=dt*(nu*Lap(V)-Adv(V,V));
+            RHS=1/dt*Div(dV);
+            std::cout<<"RHS Done"<<std::endl;
+            DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+            auto Pressure_Poisson = Lap(P);
+            auto D_y=Dy(P);
+            auto D_x=Dx(P);
+            Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
+            Solver.impose(D_y, up_p,0);
+            Solver.impose(D_x, r_p, 0);
+            Solver.impose(-D_y, dw_p,0);
+            Solver.impose(-D_x, l_p,0);
+            Solver.solve(P);
+            std::cout << "Poisson Solved" << std::endl;
+            Vtemp = V + (dV - dt*Grad(P));
+            V = Vtemp;
             divV = Div(V);
             for(int j=0;j<up_p.size();j++)
             {   auto p=up_p.get<0>(j);