diff --git a/demo/simple/1.in b/demo/simple/1.in index 7fc02d8d3cb9a9dfb28d8a40fefde6e6bddf4f29..493b429431468e07cea1f8b95e376ac5abe082a2 100644 --- a/demo/simple/1.in +++ b/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 diff --git a/demo/simple/2.in b/demo/simple/2.in index f78da94a58ba7bfe25c4c62af52f891347f3be2d..b85fa3c697740dc4a25637ede5c784bd23354bcf 100644 --- a/demo/simple/2.in +++ b/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 diff --git a/demo/simple/3.in b/demo/simple/3.in index a8d0928ae93c7866c4c03937ed48254bb86f076b..fe2bcf85fdbe7abf62119666fc12150494bba95a 100644 --- a/demo/simple/3.in +++ b/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 diff --git a/include/simple/simple.h b/include/simple/simple.h index d40a3a02f82d361c6998ae35f969233dd44d47d0..e75b930e918a11947aadaa6674b10d73a093663f 100644 --- a/include/simple/simple.h +++ b/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 diff --git a/source/simple/abstract_parameters.cpp b/source/simple/abstract_parameters.cpp index e1a8abeb4dd25d2b42841db9e565f84e3316d879..8b212a4451fdc50568971055c13e2b45ca4cb2b4 100644 --- a/source/simple/abstract_parameters.cpp +++ b/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"); diff --git a/source/simple/simple.cpp b/source/simple/simple.cpp index bb1541cc437d918472ac560fd16abaac42ad4c32..6fd2d93bf52086ece9339f18b5d8785cf4c8621c 100644 --- a/source/simple/simple.cpp +++ b/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; } } diff --git a/source/simple/solve.cpp b/source/simple/solve.cpp index f65906047a108ad984b84e34eaa1252992776f10..d68229885bc484de34f6c54c33fb5a52e49a7863 100644 --- a/source/simple/solve.cpp +++ b/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,26 +78,26 @@ 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 if (max_residual < _parameters->p_max_residual) { - return max_e; + return max_e; } 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,24 +185,25 @@ 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; while(true) { 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; }