Skip to content
Snippets Groups Projects
Commit 26f098cd authored by Kyrylo Sovailo's avatar Kyrylo Sovailo :crossed_swords:
Browse files

Fixed logics?

parent 58906b79
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
include/simple/simple.h
+ 11
− 7
View file @ 26f098cd
...@@ -56,8 +56,9 @@ namespace simple ...@@ -56,8 +56,9 @@ namespace simple
std::vector<FaceConstants> face_constants; ///< Precalculated geometry data std::vector<FaceConstants> face_constants; ///< Precalculated geometry data
//Variables //Variables
double velocity; ///< Estimate of velocity v* double velocity; ///< Current most latest estimate of velocity V
double new_velocity; ///< New estimate of velocity v*, used in Jacobi algorithm double velocity_estimate; ///< Current estimate of velocity V*
double new_velocity; ///< New estimate of velocity V*, used in Jacobi algorithm
//Cache //Cache
double a0_tilde; ///< a0_tilde = sum(a) coefficient, calculated during V* solver to be used in P' solver double a0_tilde; ///< a0_tilde = sum(a) coefficient, calculated during V* solver to be used in P' solver
...@@ -88,7 +89,8 @@ namespace simple ...@@ -88,7 +89,8 @@ namespace simple
std::vector<FaceConstants> face_constants; ///< Precalculated for each face std::vector<FaceConstants> face_constants; ///< Precalculated for each face
//Variables //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 pressure_correction; ///< Pressure correction P'
double new_pressure_correction; ///< New pressure correction P', used in Jacobi algorithm 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 double e; ///< e = sum(rho * L * v) coefficient, precalculated for each P' solution
...@@ -192,10 +194,12 @@ namespace simple ...@@ -192,10 +194,12 @@ namespace simple
double _a(double pe) const; ///< Calculates A(Pe) double _a(double pe) const; ///< Calculates A(Pe)
void _setup_fixed_pressure(); ///< Setup pressure on all cells void _setup_fixed_pressure(); ///< Setup pressure on all cells
void _setup_fixed_velocity(); ///< Setup velocity on all fixed faces void _setup_fixed_velocity(); ///< Setup velocity on all fixed faces
double _solve_pressure_correction(); ///< Solves linear system for P' void _solve_velocity(); ///< Solves linear system for V
void _solve_velocity_estimate(); ///< Solves linear system for V*, returns e double _solve_pressure_correction(); ///< Solves linear system for P', returns e
void _correct_pressure(); ///< Updates P' void _set_correct_pressure(); ///< Updates P = P* + P'
void _correct_velocity(); ///< Updates V* void _set_correct_velocity(); ///< Updates V = V* + V'
void _set_pressure_estimate(); ///< Updates V* = V
void _set_velocity_estimate(); ///< Updates P* = P
public: public:
Solver(const AbstractParameters *parameters, std::string filename); ///< Creates solver Solver(const AbstractParameters *parameters, std::string filename); ///< Creates solver
......
source/simple/solve.cpp
+ 158
− 101
View file @ 26f098cd
...@@ -23,7 +23,7 @@ void simple::Solver::_setup_fixed_pressure() ...@@ -23,7 +23,7 @@ void simple::Solver::_setup_fixed_pressure()
for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++) for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{ {
Cell *cell = *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() ...@@ -32,7 +32,98 @@ void simple::Solver::_setup_fixed_velocity()
for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++) for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{ {
Face *face = *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() ...@@ -84,25 +175,9 @@ double simple::Solver::_solve_pressure_correction()
for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++) for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{ {
Cell *cell = *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; //Do nothing at fixed pressure or fixed e
perimeter += cell->faces[n]->length;
}
}
(_parameters->p_gauss ? cell->pressure_correction : cell->new_pressure_correction) = sum / perimeter;
} }
else else
{ {
...@@ -143,114 +218,92 @@ double simple::Solver::_solve_pressure_correction() ...@@ -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(); for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
unsigned int fail_count = 0;
unsigned int iteration = 0;
while(true)
{ {
iteration++; Cell *cell = *icell;
double max_residual = -std::numeric_limits<double>::infinity(); if (cell->fixed_pressure || cell->fixed_e)
//Main calculation
for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{ {
Face *face = *iface; // Do nothing
if (face->fixed_velocity) }
else
{ {
//Do nothing at fixed velocity // Normal update
cell->pressure = cell->pressure_estimate + _parameters->beta * cell->pressure_correction;
} }
else if (face->cells.size() < 2) }
for (std::set<Cell*>::iterator icell = _cells.begin(); icell != _cells.end(); icell++)
{ {
//If it's not fixed velocity, but an edge, just copy the velocity from the neighboring cell Cell *cell = *icell;
//Because edge cells don't have area, and don't have pressure difference if (!cell->fixed_pressure && cell->fixed_e)
for (unsigned int n = 0; n < face->neighbors.size(); n++)
{ {
if (face->face_constants[n].neighbor_velocity_inwards_sign != 0) //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++)
{ {
(_parameters->v_gauss ? face->velocity : face->new_velocity) = face->face_constants[n].neighbor_velocity_sign * face->neighbors[n]->velocity; if (cell->neighbors[n] != nullptr)
} {
sum += cell->neighbors[n]->pressure_correction * cell->faces[n]->length;
perimeter += cell->faces[n]->length;
} }
} }
else cell->pressure = sum / perimeter;
{
//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) void simple::Solver::_set_correct_velocity()
if (!_parameters->v_gauss)
{ {
for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++) for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
{ {
Face *face = *iface; Face *face = *iface;
if (face->fixed_velocity) continue; //Skip fixed cells if (face->fixed_velocity || face->neighbors.size() < 2)
face->velocity = face->new_velocity; {
// Do nothing
}
else
{
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;
} }
} }
//Logics for (std::set<Face*>::iterator iface = _faces.begin(); iface != _faces.end(); iface++)
if (max_residual < _parameters->v_max_residual)
{ {
break; Face *face = *iface;
} if (!face->fixed_velocity && face->neighbors.size() < 2)
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"); // Velocity on edges is a copy of neighbor velocity
else if (_parameters->v_reaction == Reaction::warning) std::cout << "V* solver exceeded number of iterations\n"; for (unsigned int n = 0; n < face->neighbors.size(); 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"); 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++) for (std::set<Cell*>::iterator i = _cells.begin(); i != _cells.end(); i++)
{ {
Cell *cell = *i; Cell *cell = *i;
if (cell->fixed_pressure) continue; 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++) for (std::set<Face*>::iterator i = _faces.begin(); i != _faces.end(); i++)
{ {
Face *face = *i; Face *face = *i;
if (face->fixed_velocity || face->neighbors.size() < 2) continue; //Skip fixed velocity cells and edge cells if (face->fixed_velocity) 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_estimate = face->velocity;
} }
} }
...@@ -266,10 +319,9 @@ void simple::Solver::step() ...@@ -266,10 +319,9 @@ void simple::Solver::step()
while(true) while(true)
{ {
iteration++; iteration++;
_solve_velocity_estimate(); double max_e = _solve_pressure_correction(); //Solve P' (based on V)
double max_e = _solve_pressure_correction(); _set_correct_pressure(); //P = P* + P'
_correct_pressure(); _set_correct_velocity(); //V = V* + V'
_correct_velocity();
if (max_e < _parameters->max_e) if (max_e < _parameters->max_e)
{ {
std::streamsize precision = std::cout.precision(4); std::streamsize precision = std::cout.precision(4);
...@@ -287,6 +339,11 @@ void simple::Solver::step() ...@@ -287,6 +339,11 @@ void simple::Solver::step()
previous_max_e = max_e; previous_max_e = max_e;
} }
//Solution
_solve_velocity(); //Solve V (based on P)
_set_pressure_estimate(); //P* = P
_set_velocity_estimate(); //V' = V
//Progressing //Progressing
_time += _parameters->step; _time += _parameters->step;
} }
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment