diff --git a/include/simple/simple.h b/include/simple/simple.h
index b24fc5240e58fdfa59a513de50c576d028127db5..9e6955098db88ec7adc804869ee8dac3b00d0282 100644
--- a/include/simple/simple.h
+++ b/include/simple/simple.h
@@ -56,8 +56,9 @@ namespace simple
         std::vector<FaceConstants> face_constants;   ///< Precalculated geometry data
 
         //Variables
-        double velocity;                ///< Estimate of velocity v*
-        double new_velocity;            ///< New estimate of velocity v*, used in Jacobi algorithm
+        double velocity;                ///< Current most latest estimate of velocity V
+        double velocity_estimate;       ///< Current estimate of velocity V*
+        double new_velocity;            ///< New estimate of velocity V*, used in Jacobi algorithm
 
         //Cache
         double a0_tilde;                ///< a0_tilde = sum(a) coefficient, calculated during V* solver to be used in P' solver
@@ -88,7 +89,8 @@ namespace simple
         std::vector<FaceConstants> face_constants;  ///< Precalculated for each face
 
         //Variables
-        double pressure;                ///< Current estimate of pressure P*
+        double pressure;                ///< Current most latest estimate of pressure P
+        double pressure_estimate;       ///< Current estimate of pressure P*
         double pressure_correction;     ///< Pressure correction P'
         double new_pressure_correction; ///< New pressure correction P', used in Jacobi algorithm
         double e;                       ///< e = sum(rho * L * v) coefficient, precalculated for each P' solution
@@ -192,10 +194,12 @@ namespace simple
         double _a(double pe) const;                                                                                 ///< Calculates A(Pe)
         void _setup_fixed_pressure();                                                                               ///< Setup pressure on all cells
         void _setup_fixed_velocity();                                                                               ///< Setup velocity on all fixed faces
-        double _solve_pressure_correction();                                                                        ///< Solves linear system for P'
-        void _solve_velocity_estimate();                                                                            ///< Solves linear system for V*, returns e
-        void _correct_pressure();                                                                                   ///< Updates P'
-        void _correct_velocity();                                                                                   ///< Updates V*
+        void _solve_velocity();                                                                                     ///< Solves linear system for V
+        double _solve_pressure_correction();                                                                        ///< Solves linear system for P', returns e
+        void _set_correct_pressure();                                                                               ///< Updates P = P* + P'
+        void _set_correct_velocity();                                                                               ///< Updates V = V* + V'
+        void _set_pressure_estimate();                                                                              ///< Updates V* = V
+        void _set_velocity_estimate();                                                                              ///< Updates P* = P
 
     public:
         Solver(const AbstractParameters *parameters, std::string filename); ///< Creates solver
diff --git a/source/simple/solve.cpp b/source/simple/solve.cpp
index b400f5676be5b6e4e806491834fb2e1c8dc877ec..fcb71db87921d82192efe69079cb4271a39c5800 100644
--- a/source/simple/solve.cpp
+++ b/source/simple/solve.cpp
@@ -23,7 +23,7 @@ void simple::Solver::_setup_fixed_pressure()
     for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
     {
         Cell *cell = *icell;
-        if (cell->fixed_pressure) cell->pressure = _parameters->pressure(cell->center, _time);
+        if (cell->fixed_pressure) cell->pressure = cell->pressure_estimate = _parameters->pressure(cell->center, _time);
     }
 }
 
@@ -32,7 +32,98 @@ void simple::Solver::_setup_fixed_velocity()
     for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
     {
         Face *face = *iface;
-        if (face->fixed_velocity) face->velocity = _parameters->velocity(face->center, _time).dot(face->normal);
+        if (face->fixed_velocity) face->velocity = face->velocity_estimate = _parameters->velocity(face->center, _time).dot(face->normal);
+    }
+}
+
+void simple::Solver::_solve_velocity()
+{
+    double previous_residual = std::numeric_limits<double>::infinity();
+    unsigned int fail_count = 0;
+    unsigned int iteration = 0;
+    while(true)
+    {
+        iteration++;
+        double max_residual = -std::numeric_limits<double>::infinity();
+        
+        //Main calculation
+        for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+        {
+            Face *face = *iface;
+            if (face->fixed_velocity)
+            {
+                //Do nothing at fixed velocity
+            }
+            else if (face->cells.size() < 2)
+            {
+                //If it's not fixed velocity, but an edge, just copy the velocity from the neighboring cell
+                //Because edge cells don't have area, and don't have pressure difference
+                for (unsigned int n = 0; n < face->neighbors.size(); n++)
+                {
+                    if (face->face_constants[n].neighbor_velocity_inwards_sign != 0)
+                    {
+                        (_parameters->v_gauss ? face->velocity : face->new_velocity) = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
+                    }
+                }
+            }
+            else
+            {
+                //Average cell, normal calculation with friction
+                double a0 = _parameters->density * face->area / _parameters->step;  //a0 becomes a0_tilde till the end of the iteration
+                double b = a0 * face->velocity;                                     //b becomes a0_tilde*next_velocity till the end of the iteration
+                b -= face->pressure_sign * (face->cells[1]->pressure - face->cells[0]->pressure) * face->area / face->cell_cell_distance;
+                for (unsigned int n = 0; n < face->neighbors.size(); n++)
+                {
+                    //Precalculating
+                    //face velocity in the right direction, used as scalar
+                    const double neighbor_velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
+                    //face velocity in inwards direction, used as velocity
+                    const double neighbor_velocity_inwards = face->face_constants[n].neighbor_velocity_inwards_sign * neighbor_velocity;
+                    const double face_face_distance = face->face_constants[n].face_face_distance;
+                    const double border_length = face->face_constants[n].border_length;
+
+                    //Kernel
+                    const double d = _parameters->viscosity / face_face_distance;
+                    const double pe = _parameters->density * neighbor_velocity_inwards * face_face_distance / _parameters->viscosity;
+                    const double f = neighbor_velocity_inwards * _parameters->density;
+                    const double a = d * border_length * _a(pe) + std::max(0.0, f * border_length);
+                    a0 += a;
+                    b += a * neighbor_velocity;
+                }
+                face->a0_tilde = a0;
+                const double residual = std::abs(face->velocity * a0 - b);
+                if (residual > max_residual) max_residual = residual;
+                (_parameters->v_gauss ? face->velocity : face->new_velocity) = b / a0;
+            }
+        }
+
+        //Applying new_velocity (Jacobi algorithm)
+        if (!_parameters->v_gauss)
+        {
+            for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+            {
+                Face *face = *iface;
+                if (face->fixed_velocity) continue; //Skip fixed cells
+                face->velocity = face->new_velocity;
+            }
+        }
+
+        //Logics
+        if (max_residual < _parameters->v_max_residual)
+        {
+            break;
+        }
+        if (iteration > _parameters->v_max_iteration)
+        {
+            if (_parameters->v_reaction == Reaction::error) throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Number of iterations exceeded");
+            else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations\n";
+            break;
+        }
+        if (max_residual > previous_residual)
+        {
+            if (fail_count++ > 100) break; //throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Residual does not converge");
+        }
+        previous_residual = max_residual;
     }
 }
 
@@ -84,25 +175,9 @@ double simple::Solver::_solve_pressure_correction()
         for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
         {
             Cell *cell = *icell;
-            if (cell->fixed_pressure)
+            if (cell->fixed_pressure || cell->fixed_e)
             {
-                //Do nothing at fixed pressure
-            }
-            else if (cell->fixed_e)
-            {
-                //All neighboring cells are fixed velocity and/or edge
-                //Pressure of the cell is a sum of neighboring cells
-                double sum = 0;
-                double perimeter = 0;
-                for (unsigned int n = 0; n < cell->neighbors.size(); n++)
-                {
-                    if (cell->neighbors[n] != nullptr)
-                    {
-                        sum += cell->neighbors[n]->pressure_correction * cell->faces[n]->length;
-                        perimeter += cell->faces[n]->length;
-                    }
-                }
-                (_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) = sum / perimeter;
+                //Do nothing at fixed pressure or fixed e
             }
             else
             {
@@ -143,114 +218,92 @@ double simple::Solver::_solve_pressure_correction()
     }
 }
 
-void simple::Solver::_solve_velocity_estimate()
+void simple::Solver::_set_correct_pressure()
 {
-    double previous_residual = std::numeric_limits<double>::infinity();
-    unsigned int fail_count = 0;
-    unsigned int iteration = 0;
-    while(true)
+    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
     {
-        iteration++;
-        double max_residual = -std::numeric_limits<double>::infinity();
-        
-        //Main calculation
-        for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+        Cell *cell = *icell;
+        if (cell->fixed_pressure || cell->fixed_e)
         {
-            Face *face = *iface;
-            if (face->fixed_velocity)
-            {
-                //Do nothing at fixed velocity
-            }
-            else if (face->cells.size() < 2)
-            {
-                //If it's not fixed velocity, but an edge, just copy the velocity from the neighboring cell
-                //Because edge cells don't have area, and don't have pressure difference
-                for (unsigned int n = 0; n < face->neighbors.size(); n++)
-                {
-                    if (face->face_constants[n].neighbor_velocity_inwards_sign != 0)
-                    {
-                        (_parameters->v_gauss ? face->velocity : face->new_velocity) = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
-                    }
-                }
-            }
-            else
-            {
-                //Average cell, normal calculation with friction
-                double a0 = _parameters->density * face->area / _parameters->step;  //a0 becomes a0_tilde till the end of the iteration
-                double b = a0 * face->velocity;                                     //b becomes a0_tilde*next_velocity till the end of the iteration
-                b -= face->pressure_sign * (face->cells[1]->pressure - face->cells[0]->pressure) * face->area / face->cell_cell_distance;
-                for (unsigned int n = 0; n < face->neighbors.size(); n++)
-                {
-                    //Precalculating
-                    //face velocity in the right direction, used as scalar
-                    const double neighbor_velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
-                    //face velocity in inwards direction, used as velocity
-                    const double neighbor_velocity_inwards = face->face_constants[n].neighbor_velocity_inwards_sign * neighbor_velocity;
-                    const double face_face_distance = face->face_constants[n].face_face_distance;
-                    const double border_length = face->face_constants[n].border_length;
-
-                    //Kernel
-                    const double d = _parameters->viscosity / face_face_distance;
-                    const double pe = _parameters->density * neighbor_velocity_inwards * face_face_distance / _parameters->viscosity;
-                    const double f = neighbor_velocity_inwards * _parameters->density;
-                    const double a = d * border_length * _a(pe) + std::max(0.0, f * border_length);
-                    a0 += a;
-                    b += a * neighbor_velocity;
-                }
-                face->a0_tilde = a0;
-                const double residual = std::abs(face->velocity * a0 - b);
-                if (residual > max_residual) max_residual = residual;
-                (_parameters->v_gauss ? face->velocity : face->new_velocity) = b / a0;
-            }
+            // Do nothing
         }
+        else
+        {
+            // Normal update
+            cell->pressure = cell->pressure_estimate + _parameters->beta * cell->pressure_correction;
+        }
+    }
 
-        //Applying new_velocity (Jacobi algorithm)
-        if (!_parameters->v_gauss)
+    for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
+    {
+        Cell *cell = *icell;
+        if (!cell->fixed_pressure && cell->fixed_e)
         {
-            for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+            //Pressure of the cell is a sum of neighboring cells
+            double sum = 0;
+            double perimeter = 0;
+            for (unsigned int n = 0; n < cell->neighbors.size(); n++)
             {
-                Face *face = *iface;
-                if (face->fixed_velocity) continue; //Skip fixed cells
-                face->velocity = face->new_velocity;
+                if (cell->neighbors[n] != nullptr)
+                {
+                    sum += cell->neighbors[n]->pressure_correction * cell->faces[n]->length;
+                    perimeter += cell->faces[n]->length;
+                }
             }
+            cell->pressure = sum / perimeter;
         }
+    }
+}
 
-        //Logics
-        if (max_residual < _parameters->v_max_residual)
+void simple::Solver::_set_correct_velocity()
+{
+    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+    {
+        Face *face = *iface;
+        if (face->fixed_velocity || face->neighbors.size() < 2)
         {
-            break;
+            // Do nothing
         }
-        if (iteration > _parameters->v_max_iteration)
+        else
         {
-            if (_parameters->v_reaction == Reaction::error) throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Number of iterations exceeded");
-            else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations\n";
-            break;
+            face->velocity = face->velocity_estimate - face->pressure_sign * (face->cells[1]->pressure_correction - face->cells[0]->pressure_correction) * face->area / face->cell_cell_distance / face->a0_tilde;
         }
-        if (max_residual > previous_residual)
+    }
+
+    for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
+    {
+        Face *face = *iface;
+        if (!face->fixed_velocity && face->neighbors.size() < 2)
         {
-            if (fail_count++ > 100) break; //throw std::runtime_error("simple::Solver::_solve_velocity_estimate(): Residual does not converge");
+            // Velocity on edges is a copy of neighbor velocity
+            for (unsigned int n = 0; n < face->neighbors.size(); n++)
+            {
+                if (face->face_constants[n].neighbor_velocity_inwards_sign != 0)
+                {
+                    face->velocity = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity;
+                }
+            }
         }
-        previous_residual = max_residual;
     }
 }
 
-void simple::Solver::_correct_pressure()
+void simple::Solver::_set_pressure_estimate()
 {
     for (std::set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
     {
         Cell *cell = *i;
         if (cell->fixed_pressure) continue;
-        cell->pressure = cell->pressure + _parameters->beta * cell->pressure_correction;
+        cell->pressure_estimate = cell->pressure;
     }
 }
 
-void simple::Solver::_correct_velocity()
+void simple::Solver::_set_velocity_estimate()
 {
     for (std::set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
     {
         Face *face = *i;
-        if (face->fixed_velocity || face->neighbors.size() < 2) continue; //Skip fixed velocity cells and edge cells
-        face->velocity = face->velocity - face->pressure_sign * (face->cells[1]->pressure_correction - face->cells[0]->pressure_correction) * face->area / face->cell_cell_distance / face->a0_tilde;
+        if (face->fixed_velocity) continue;
+        face->velocity_estimate = face->velocity;
     }
 }
 
@@ -266,10 +319,9 @@ void simple::Solver::step()
     while(true)
     {
         iteration++;
-        _solve_velocity_estimate();
-        double max_e = _solve_pressure_correction();
-        _correct_pressure();
-        _correct_velocity();
+        double max_e = _solve_pressure_correction(); //Solve P' (based on V)
+        _set_correct_pressure(); //P = P* + P'
+        _set_correct_velocity(); //V = V* + V'
         if (max_e < _parameters->max_e)
         {
             std::streamsize precision = std::cout.precision(4);
@@ -287,6 +339,11 @@ void simple::Solver::step()
         previous_max_e = max_e;
     }
 
+    //Solution
+    _solve_velocity(); //Solve V (based on P)
+    _set_pressure_estimate(); //P* = P
+    _set_velocity_estimate(); //V' = V
+
     //Progressing
     _time += _parameters->step;
 }
\ No newline at end of file