Potential bugfixes

demo/simple/1.in

@@ -10,17 +10,17 @@ METHOD LINEAR
 DENSITY 1000
 VISCOSITY 1
 STEP 0.01
-BETA 1
+BETA 0.001
 MAX_ITERATION 1000
-MAX_B 0.001
+MAX_E 0.001
 P_MAX_RESIDUAL 0
-P_MAX_ITERATION 5
+P_MAX_ITERATION 20000
 P_SOLVER GAUSS
-P_REACTION WARNING
+P_REACTION IGNORE
 V_MAX_RESIDUAL 0
-V_MAX_ITERATION 5
+V_MAX_ITERATION 1000
 V_SOLVER GAUSS
-V_REACTION WARNING
+V_REACTION IGNORE
 
 STEPS 100
 DIVIDER 50
\ No newline at end of file

demo/simple/2.in

@@ -10,7 +10,7 @@ METHOD LINEAR
 DENSITY 1000
 VISCOSITY 1
 STEP 0.01
-BETA 1
+BETA 0.5
 MAX_ITERATION 1000
 MAX_B 0.001
 P_MAX_RESIDUAL 0

demo/simple/3.in

@@ -9,7 +9,7 @@ METHOD LINEAR
 DENSITY 1000
 VISCOSITY 1
 STEP 0.01
-BETA 1
+BETA 0.5
 MAX_ITERATION 1000
 MAX_B 0.001
 P_MAX_RESIDUAL 0

include/simple/simple.h

@@ -140,8 +140,8 @@ namespace simple
         bool exact;                     ///< Potential function is valid on entire field, used for testing
         //Iterative solver parameters
         double beta;                    ///< Undershooting/overshooting parameter
-        double max_iteration;           ///< Maximum iteration
-        unsigned int max_b;             ///< Maximum b (flow in/out of cells)
+        double max_e;                   ///< Maximum e (flow in/out of cells)
+        unsigned int  max_iteration;    ///< Maximum iteration
         //Linear solver for P' parameters
         double p_max_residual;          ///< Maximum residual for P' solver
         unsigned int p_max_iteration;   ///< Maximum iteration for P' solver

source/simple/abstract_parameters.cpp

@@ -24,7 +24,7 @@ simple::AbstractParameters::AbstractParameters(numsimulation::ParameterReader &r
     //Iterative solver parameters
     beta = reader.get_real("BETA");
     max_iteration = reader.get_real("MAX_ITERATION");
-    max_b = reader.get_real("MAX_B");
+    max_e = reader.get_real("MAX_E");
     //Linear solver for P' parameters
     p_max_residual = reader.get_real("P_MAX_RESIDUAL");
     p_max_iteration = reader.get_integer("P_MAX_ITERATION");

source/simple/simple.cpp

@@ -513,6 +513,7 @@ void simple::Solver::_setup_pressure()
     {
         Cell *cell = *icell;
         cell->pressure = _parameters->pressure(cell->center, 0);
+        cell->pressure_correction = 0;
     }
 }
 

source/simple/solve.cpp

@@ -38,7 +38,7 @@ void simple::Solver::_setup_fixed_velocity()
 
 double simple::Solver::_solve_pressure_correction()
 {
-    //Precalculating on faces
+    //Precalculating c on faces
     for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
     {
         Face *face = *iface;
@@ -46,12 +46,15 @@ double simple::Solver::_solve_pressure_correction()
         face->c = face->area / (face->a0_tilde * face->length);
     }
 
-    //Precalculating on cells
+    //Precalculating e and c0_tilde on cells, setting pressure correction to zero
     double max_e = -std::numeric_limits<double>::infinity();
     for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
     {
         Cell *cell = *icell;
         if (cell->fixed_pressure) continue; //Skip fixed cells
+        
+        cell->pressure_correction = 0;
+
         cell->e = 0;
         cell->c0_tilde = 0;
         for (unsigned int f = 0; f < cell->faces.size(); f++)
@@ -75,15 +78,15 @@ double simple::Solver::_solve_pressure_correction()
             Cell *cell = *icell;
             if (cell->fixed_pressure) continue; //Skip fixed cells
 
-            double sum = cell->e;
+            double sum = 0;
             for (unsigned int f = 0; f < cell->faces.size(); f++)
             {
                 if (cell->neighbors[f] == nullptr) continue; //Skip if no cell on other side
                 sum += cell->faces[f]->c * cell->neighbors[f]->pressure_correction;
             }
-            double residual = std::abs(sum - cell->c0_tilde * cell->pressure_correction);
+            double residual = std::abs(cell->e + sum - cell->c0_tilde * cell->pressure_correction);
             if (residual > max_residual) max_residual = residual;
-            (_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) = sum / cell->c0_tilde;
+            (_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) = (cell->e + sum) / cell->c0_tilde;
         }
 
         //Logics
@@ -94,7 +97,7 @@ double simple::Solver::_solve_pressure_correction()
         if (iteration > _parameters->p_max_iteration)
         {
             if (_parameters->p_reaction == Reaction::error) throw std::runtime_error("void simple::Solver::_solve_pressure_correction(): Number of iterations exceeded");
-            else if (_parameters->p_reaction == Reaction::warning) std::cout << "P' solver exceeded number of iterations";
+            else if (_parameters->p_reaction == Reaction::warning) std::cout << "P' solver exceeded number of iterations\n";
             return max_e;
         }
     }
@@ -160,7 +163,7 @@ void simple::Solver::_solve_velocity_estimate()
         if (iteration > _parameters->v_max_iteration)
         {
             if (_parameters->v_reaction == Reaction::error) throw std::runtime_error("void simple::Solver::_solve_velocity_estimate(): Number of iterations exceeded");
-            else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations";
+            else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations\n";
             break;
         }
     }
@@ -182,13 +185,14 @@ void simple::Solver::_correct_velocity()
     {
         Face *face = *i;
         if (face->cells.size() < 2) continue;
-        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;
+        face->velocity = face->velocity + _parameters->beta * face->pressure_sign * (face->cells[1]->pressure_correction - face->cells[0]->pressure_correction) * face->area / face->cell_cell_distance / face->a0_tilde;
     }
 }
 
 void simple::Solver::step()
 {
-
+    _setup_fixed_pressure();
+    _setup_fixed_velocity();
 
     //Iteration
     unsigned int iteration = 0;
@@ -196,10 +200,10 @@ void simple::Solver::step()
     {
         iteration++;
         _solve_velocity_estimate();
-        double max_b = _solve_pressure_correction();
+        double max_e = _solve_pressure_correction();
         _correct_pressure();
         _correct_velocity();
-        if (max_b < _parameters->max_b)
+        if (max_e < _parameters->max_e)
         {
             break;
         }